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Preface 


The  study  documented  herein  was  undertaken  as  part  of  the  Work  Unit 
“Innovative  Geophysical  Technologies  for  Enhanced  Buried  Unexploded 
Ordnance  (UXO)  Discrimination”  (AF25,  6.2)  and  funded  by  the  U.S.  Army 
Engineer  Research  and  Development  Center  (ERDC),  Vicksburg,  MS,  under  the 
Environmental  Quality  Technology  (EQT)  Program.  Dr.  Dwain  K.  Butler, 
Geotechnical  and  Structures  Laboratory  (GSL),  is  Principal  Investigator  for  the 
work  unit,  Dr.  John  M.  Cullinane,  Environmental  Laboratory  (EL),  is  EQT 
Program  Manager,  and  Dr.  Ernesto  R.  Cespedes,  EL,  is  UXO  Focus  Area 
Manager  for  EQT.  Dr.  David  W.  Pittman  was  Acting  Director,  GSL.  The  work 
was  performed  at  the  University  of  British  Columbia  (UBC),  Vancouver,  British 
Columbia,  Canada.  Funding  of  the  UBC  effort  was  by  means  of  Army  Research 
Office  Grant  41 262-EV.  Dr.  Russell  Harmon,  Terrestrial  Sciences  Program 
Manager,  administered  the  UBC  grant. 

Dr.  Stephen  D.  Billings,  Mr.  Leonard  R.  Pasion,  and  Dr.  Douglas  W. 
Oldenburg,  UBC,  performed  the  work  during  the  period  1  April  2001  through 
3 1  March  2002.  The  work  exploits  a  forward  model,  based  on  the  multipole 
expansion  of  the  induced  magnetic  field  in  a  spheroid  model  of  UXO,  in  an 
inversion  algorithm  to  recover  the  parameters  of  a  spheroid  model  fit  to  measured 
total  field  magnetic  data.  Prolate  spheroids  are  proposed  as  an  approximate 
model  for  UXO  by  a  several  authors  (McFee  1989;  Altshuler  1996;  Butler  et  al. 
1998).  The  total  magnetic  field  inversion  algorithm  is  applied  to  datasets  from 
Guthrie  Road,  Montana,  and  Former  Fort  Ord,  California.  Based  on  results  on 
the  inversions,  discrimination  and  identification  approaches  are  proposed. 
Application  of  the  discrimination  and  identification  algorithms  to  total  field 
magnetic  data  illustrates  the  potential  to  significantly  reduce  the  cost  of  UXO 
cleanup  by  discrimination  of  false  alarm  anomalies  from  UXO  anomalies. 

At  the  time  of  publication  of  this  report,  Dr.  James  R.  Houston  was  Director, 
ERDC.  Commander  and  Executive  Director  of  ERDC  was  COL  John  W. 

Morris  III,  EN. 


The  contents  of  this  report  are  not  to  be  used  for  advertising,  publication,  or 
promotional  purposes.  Citation  of  trade  names  does  not  constitute  an  official 
endorsement  or  approval  of  the  use  of  such  commercial  products. 


1  Introduction 


Unexploded  ordnance  pose  a  significant  public  safety  hazard  in  many  parts  of  the 
world.  They  occur  on  or  near  the  surface  and  down  to  depths  of  several  meters  below  the 
ground.  One  source  of  UXO’s  are  armed  conflicts  both  old  and  more  recent  in  regions 
such  as  Europe  (World  Wars  I  and  II),  South-East  Asia,  Afghanistan,  Yugoslavia  and 
parts  of  Africa.  They  are  also  a  significant  problem  in  countries  such  as  Canada  and 
the  United  States,  where  they  are  present  in  areas  used  for  military  training  and  firing 
ranges;  largely  as  vestiges  of  World  War  II  and  the  Cold  War.  It  is  estimated  that  15 
million  acres  of  the  United  States  are  contaminated  by  UXO  (FAC,  1996),  with  a  clean¬ 
up  time-frame  in  the  decades  and  a  staggering  cost  estimate,  using  existing  technology, 
in  the  tens  to  hundreds  of  billions  of  dollars  (GAO,  2001). 

The  most  well  established  techniques  for  ordnance  detection  are  magnetics  and  elec¬ 
tromagnetics  (Bulter  et  al.,  1998).  These  methods  are  very  effective  at  locating  buried 
metallic  objects  such  as  UXO.  However,  discriminating  between  intact  UXOs  and  unin¬ 
teresting  objects  such  as  metallic  debris  and  shrapnel  is  a  significant  challenge  (Bulter  et 
al.,  1998).  Often  little  or  no  effort  is  invested  in  discrimination  and,  consequently,  many 
holes  are  excavated  for  each  ordnance  item  recovered.  For  example,  Putnam  (2001)  re¬ 
ports  that  from  49,521  anomalies  excavated  at  Kaho’olawe,  Hawaii,  only  3%  (or  1,486) 
turned  out  to  be  UXO.  In  this  report  we  analyze  magnetics  data  collected  at  Guthrie 
Road,  Helana  Valley,  Montana  (Youmans  and  Daehn,  1999).  840  holes  were  excavated 
with  83  UXOs  recovered,  thus  there  was  a  10%  success  rate.  It  is  perhaps  a  reflection 
of  the  state  of  our  ability  to  discriminate  that  this  survey  was  considered  a  success  by 
all  parties  involved. 

The  simplest  form  of  discrimination  is  to  threshold  the  collected  data  in  some  way 
(for  magnetics  this  is  usually  achieved  by  first  calculating  the  analytic  signal  and  then 
thresholding,  Roest  et  al.  (2001)).  As  the  threshold  level  is  decreased,  the  number  of 
items  that  need  to  be  excavated  increases.  As  demonstrated  by  the  Ordnance  Detection 
and  Discrimination  Study  at  the  Former  Fort  Ord  (Asch  and  Staes,  2001),  for  such 
techniques  to  yield  a  high  probability  of  detection,  a  very  high  false  alarm  rate  results. 

Location  of  UXOs  involves  the  successful  detection  and  identification  of  compact 
metallic  objects.  Therefore,  a  better  method  for  discrimination  is  to  parameterize  the 
characteristics  of  the  target  in  some  way  (e.g.  by  its  location,  orientation,  shape,  size 
and  material  properties)  and  develop  a  good  forward  model  that  allows  predicted  data 
to  be  generated  from  a  given  target  realization.  Inversion  can  then  be  used  to  find  the 
set  of  target  parameters  that  best  fits  the  data.  The  target  parameters  are  then  used 
to  make  decisions  regarding  discrimination  (UXO  versus  non-UXO)  and  identification 
(what  is  the  ordnance  type). 

A  previous  report  (Pasion  and  Oldenburg,  2001)  considered  the  development  of  an 
inversion  methodology  for  time-domain  electromagnetic  data.  In  this  report,  we  con¬ 
sider  the  development  of  a  methodology  for  total-field  magnetics  data.  Previous  work 
in  this  area  has  been  conducted  by  (amongst  others)  the  Canadian  Defence  Research 
Establishment  (McFee,  1989)  and  the  Naval  Research  Laboratory  (Nelson  et  al.,  1998) 
and  private  companies. 

Spheroids  have  been  proposed  as  an  approximate  parameterization  of  ordnance  by 
several  authors  (McFee,  1989;  Altshuler,  1996;  Bulter  et  al.,  1998).  While  the  spheroid 
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does  not  capture  the  top-bottom  asymmetry  of  many  ordnance  items,  close  agreement 
between  observed  anomalies  over  test  stands  and  spheroid  fits  have  been  demonstrated 
(McFee,  1989;  Bulter  et  al.,  1998).  Furthermore,  the  magnetic  anomaly  from  a  solid 
spheroid  has  been  shown  to  be  very  similar  to  a  hollow  spheroid  (Altshuler,  1996). 
Therefore,  we  will  use  solid  spheroids  as  the  basis  of  our  magnetic  modelling  of  TJXO’s. 

The  response  of  a  spheroid  (in  fact  any  compact  body)  can  be  decomposed  into  a 
series  of  moments  by  a  multipole  expansion  (Stratton,  1941).  The  response  of  the  dipole 
component  dominates  the  anomalies  caused  by  most  buried  ordnance  due  to  the  rapid 
falloff  with  distance  of  the  other  components.  Consequently,  the  main  inversion  routine 
developed  in  this  report  attempts  to  recover  the  dipole  moment  and  position  that  best 
matches  an  observed  anomaly.  The  recovered  dipole  moment  is  then  used  to  rank  items 
according  to  the  likelihood  they  are  UXO.  We  prefer  such  a  ranking  over  a  definitive 
statement  such  as  ordnance  or  non-ordnance  because  it  is  not  possible  to  make  such  a 
black  and  white  distinction. 

The  recovered  dipole  moment  can  also  be  used  for  ordnance  identification.  There  is 
some  ambiguity  in  this  process  because  the  dipole  moment  is  the  product  of  magneti¬ 
zation  with  volume.  Due  to  self-demagnetization  effects  the  magnetization  can  change 
significantly  as  the  orientation  of  the  body  relative  to  the  Earth’s  field  changes.  Thus 
orientation  can  be  traded  with  volume  to  produce  the  same  or  similar  dipole  moment. 
By  recognizing  that  there  are  a  finite  number  of  ordnance  types  present  in  an  area  (we 
build  up  a  library),  this  ambiguity  can  be  reduced  significantly. 

To  attempt  to  constrain  the  orientation  we  develop  two  further  inversion  routines. 
The  first  uses  the  octupole  component  of  the  spheroid,  and  finds  the  optimum  orientation 
and  fit  for  each  of  the  items  in  our  library.  That  item  with  the  best  fit  is  assumed  to  be 
the  source  of  the  anomaly.  The  second  method  aims  to  recover  the  dipole  and  quadrupole 
moments  of  an  arbitrary  axially  symmetric  body. 

No  inversion  methodology  would  be  complete  without  an  analysis  of  the  uncertainties 
in  the  recovered  parameters.  The  last  section  of  the  report  considers  the  development  of 
a  linearized  error  analysis.  We  also  investigate  the  distribution  of  the  residuals  so  that 
we  can  test  our  assumption  of  Gaussian  statistics. 
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2  Magnetic  field  of  a  spheroid 


The  forward  modelling  task  may  be  stated  as  follows;  calculate  the  magnetic  anomaly 
at  the  arbitrary  point  (x,  y,  z)  given  a  spheroid  of  diameter  o,  aspect  ratio  e  (i.e.  length  = 
ae)  and  permeability  y  oriented  at  an  angle  of  8  degrees  to  the  horizontal  with  an 
azimuth  of  0  degrees  clockwise  from  North,  at  a  depth  z  below  the  surface  and  at 
the  horizontal  position  (x,  y).  The  geometry  of  this  situation  is  shown  in  Figure  1. 
Full  characterization  of  the  anomaly  caused  by  a  spheroid  requires  the  Earth’s  field 
plus  the  eight  parameters  ( x ,  y,  z,  8,  (/>,  a,  e,  y).  For  the  Earth’s  field  we  use  geographical 
coordinates;  Bc  =  (Box,  BOJ/,  Boz),  i.e.  x  is  positive  to  the  East,  y  is  positive  to  the 
North  and  2  is  positive  upwards.  Note  that  this  differs  from  the  definition  used  to 
specify  the  International  Geomagnetic  Reference  Field  (IGRF)  where  x  is  positive  to 
the  North,  y  is  positive  to  the  East  and  z  is  positive  downwards. 


Figure  1:  Geometry  for  modelling  the  response  of  a  spheroid. 

The  magnetic  anomaly  of  a  buried  ferrous  item  arises  from  both  remnant  M rem  and 
induced  magnetization  M,„;  the  total  magnetization  M  is  then 

M  —  Mrem  +  Mj„  (1) 

Remnant  magnetism  is  present  even  in  the  absence  of  an  inducing  field  and  is  due  to 
magnetic  moments  being  locked  into  alignment  once  the  steel  casing  cools  to  below  the 
Curie  temperature  in  the  presence  of  an  external  field.  We  will  discuss  this  type  of 
magnetism  latter  but  for  now  we  will  concentrate  on  the  induced  magnetism  that  arises 
because  magnetic  domains  in  a  ferrous  material  tend  to  align  with  the  direction  of  the 
Earth’s  field.  The  ease  with  which  the  moments  align,  and  hence  the  strength  of  the 
magnetism,  depends  on  the  magnetic  permeability  of  the  steel.  For  a  compact  body 
such  as  a  spheroid,  demagnetization  effects  become  important.  This  phenomena  refers 
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to  the  extent  that  the  induced  field  is  reduced  due  to  the  shape  of  the  spheroid.  It 
arises  due  to  the  boundary  conditions  that  the  field  must  satisfy  across  a  discontinuity 
in  magnetic  permeability.  A  solution  of  the  boundary  value  problem  (Stratton,  1941) 
shows  that  the  demagnetization  factors  are 

F ■  = _ ~  1 _  (o) 

1  1  +  -  l)/2  W 

where  /j,r  is  the  relative  permeability  (/i  =  /zrp0,  where  fi0  =  4tt  x  10-7  H/m  is  the 
permeability  of  free  space)  and  the  cq’s  are  factors  which  depend  on  the  shape  of  the 
spheroid  (see  Appendix  A).  Once  the  relative  permeability  exceeds  a  few  hundred  units, 
the  induced  magnetization  becomes  virtually  independent  of  susceptibility.  This  can 
easily  be  seen  from  Equation  2  with  large  \xr  because  then 

(3) 

This  allows  us  to  eliminate  nr  from  our  calculations  because  steel  typically  has  suscep¬ 
tibilities  of  several  hundred  to  over  a  thousand  (Lide,  2001). 

The  demagnetization  factors,  together  with  the  Earth’s  field  B0  (in  spheroid  centered 
coordinates),  determine  the  strength  of  the  induced  magnetization, 


Min  =  Mo  1FB0  (4) 

where  we  have  written  F  as  a  3  x  3  diagonal  matrix  with  (F2,  F2,  F3)  along  the  diagonal 
(remember  F\  =  F2).  An  important  consequence  of  demagnetization  is  that  the  induced 
magnetism  can  change  significantly  with  orientation.  For  instance,  with  e  =  4we  find 
F2  =  2.1  and  F3  =  12.7,  so  that  the  magnetization  when  the  spheroid  semi-major  axis 
is  aligned  with  the  field  is  approximately  6  times  greater  than  when  it  is  perpendicular. 

Equation  4  returns  Mm  in  spheroid  centered  coordinates  and  assumes  B0  is  also 
spheroid  centered.  To  use  geographical  coordinates  we  utilize  the  Euler  rotation  tensor 


A  = 


cos  ip 
sin  9  sin  ip 
cos  9  sin  ip 


—  sin^ 
sin  9  cos  ip 
cos  9  cos  ip 


0 

cos  9 
—  sin  9 


(5) 


to  rotate  B0  to  spheroid  centered  coordinates  and  then  use  its  inverse  AT  to  rotate  M,„ 
to  geographical  coordinates, 


Mjn  =  ii-1  AtFABg  (6) 

With  this  uniform  magnetization,  the  magnetic  field  of  the  spheroid  is  given  by  the 
expression 


B=rJs?M‘-dS  <7> 

where  r'  is  the  distance  from  the  source  to  the  observation  point  (Figure  1).  An  exact 
solution  of  this  equation  using  prolate  spheroidal  harmonics  is  known  (Stratton,  1941) 
but  we  chose  to  use  a  multipole  expansion;  the  details  can  be  found  in  Appendix  A  and 
McFee  (1989).  The  first  non-zero  moment  is  the  dipole  m  which  is  given  by 

m  =  VMin  =  ^ea3Min  (8) 
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where  V  is  the  volume  of  the  spheroid.  The  magnetic  field  from  the  dipole  term  is  then 

B  =  4^0  “]*-“)  W 

The  octupole  is  the  next  non-zero  moment  and  is  a  rank  3  tensor  with  27  components. 
Symmetry  reduces  the  number  of  independent  components  to  6.  The  specification  of 
the  octupole  term  is  quite  complicated,  so  we  leave  the  details  to  Appendix  A  and  just 
note  that  its  field  dies  off  as  r-5. 

The  rapid  falloff  of  the  octupole  term  with  distance  means  that  once  the  sensor 
distance  exceeds  a  few  body  lengths,  the  field  is  essentially  dipolar.  For  example,  Fig¬ 
ure  2a  plots  the  dipole  and  octupole  fields  from  a  spheroid  model  of  a  105-mm  ordnance 
(a  =  0.105  m,  e  =  3.7)  orientated  at  45°  to  the  vertical,  along  horizontal  fines  at  a  height 
of  ae  and  2ae  above  the  ordnance.  The  responses  are  normalized  by  the  maximum  dipole 
response  at  that  height.  At  a  height  of  ae  the  peak  of  the  octupole  term  is  a  little  over 
15%  of  the  peak  of  the  dipole  term.  For  the  height  2ae  the  octupole  peak  is  less  than 
5%  the  size  of  the  dipole  peak.  Figure  2b  shows  how  this  ratio  rapidly  decreases  with 
increasing  height  above  the  ordnance. 

Considering  that  real  magnetic  profiles  are  contaminated  by  noise,  observing  the 
octupole  response  in  the  presence  of  the  dipole  response  will  generally  be  very  difficult. 
It  may  be  possible  when  the  sensor  is  very  close  to  the  ordnance,  or  when  the  noise 
levels  are  low  and  the  spatial  positioning  extremely  accurate.  We  now  examine  the 
consequence  of  the  dipolar  nature  of  the  observed  fields  on  both  discrimination  and 
identification. 

Discrimination  using  magnetometry 

The  induced  dipole  moment  of  a  sphere  is  constrained  to  fie  in  the  direction  of  the 
applied  field.  On  the  other  hand,  the  angle  the  induced  dipole  moment  of  a  spheroid 
makes  with  the  Earth’s  field  will  vary  with  its  orientation.  The  angle  ip  between  the 
induced  moment  and  the  applied  field  will  be  (see  Appendix  B) 


ip  =  arccos 


F-2  sin2  9  +  F3  cos2  9 
Jf$  sin2  9  +  Ff  cos2  9 


(10) 


where  9  is  the  angle  of  the  spheroid  semi-major  axis  with  the  Earth’s  field  and  F \  and 
F3  are  the  demagnetization  factors  of  the  spheroid.  In  Figure  3a  we  show  the  variation 
in  the  angle  ip  as  the  orientation  9  changes  for  a  spheroid  with  an  aspect-ratio  of  5.  The 
angle  ip  achieves  a  maximum  value  that  depends  on  the  shape  of  the  spheroid.  For  a 
given  aspect  ratio,  the  spheroid  orientation  9  that  results  in  the  maximum  value  for  ip 
is  (see  Appendix  B) 

9  =  arctan 


Fj  -  F3F22  -  2F2Fj\ 
F2  +  F2F$-2F3F$) 


The  maximum  angle  for  a  given  aspect  ratio  is  then  obtained  by  substituting  Equa¬ 
tion  (11)  into  (10).  The  variation  of  this  maximum  angle  with  shape  is  shown  in  Fig¬ 
ure  3b.  Most  UXO’s  have  eccentricities  between  3  and  7,  so  that  55°  is  an  upper  bound 
on  the  angle.  This  suggests  a  possible  method  for  discrimination  (Altshuler,  1996). 

When  an  ordnance  hits  the  ground  but  does  not  explode,  the  impact  is  usually  suf¬ 
ficient  to  cause  shock  demagnetization;  the  shock  of  impact  causes  magnetic  domains 
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Distance  in  body  lengths 


Figure  2:  (a)  Dipole  and  octupole  terms  at  one  and  two  body  lengths  distance;  (b)  Log10 
of  the  ratio  of  the  octupole  to  dipole  terms  for  a  105-mm  projectile  at  45°  inclination 

within  the  ordnance  to  become  randomly  aligned,  thus  erasing  any  remnant  magnetiza¬ 
tion.  On  the  other  hand,  when  the  ordnance  explodes,  the  shrapnel  is  heated,  deformed 
and  cooled  in  the  presence  of  the  Earth’s  field  and  can  acquire  permanent  magnetism. 
Putting  all  of  the  above  together,  it  is  conjectured  that  ordnance  should  not  have  a 
dipole  moment  that  exceeds  an  angle  of  60°  from  the  Earth’s  field.  The  success  of  this 
method  of  discrimination  in  reducing  false  alarms  has  been  demonstrated  in  Nelson  et 
al.  (1998).  However,  they  did  note  that  subsequent  events  (such  as  lightning  strikes) 
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Figure  3:  (a)  Change  in  the  moment  angle  as  the  orientation  varies  for  an  aspect  ratio 
of  5;  and  (b)  Maximum  angle  between  the  induced  dipole  moment  of  a  spheroid  and  the 
Earth’s  field 

can  cause  the  ordnance  to  acquire  a  remnant  magnetization  and  hence  degrade  the  per¬ 
formance  of  the  discrimination  method.  Also,  the  assumption  that  ordnance  are  fully 
demagnetized  on  impact  may  be  violated  in  practice.  We  shall  return  to  this  point  later, 
but  for  now  assume  that  the  initial  conjecture  is  valid. 
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90 


1 


Figure  4:  Magnitude  of  dipole  moment  (Am2)  and  angle  from  Earth’s  field  for  60,  76, 
81,  105  and  155-mm  projectiles 


Ordnance  identification  using  magnetometry 

We  now  turn  to  the  issue  of  identifying  (or  classifying)  ordnance  items  from  their  dipole 
field.  For  a  given  spheroid,  the  induced  dipole  moment  will  be  dependent  on  the  angle  6 
that  the  spheroid  axis  makes  with  the  Earth’s  field.  There  is  no  azimuthal  dependence 
due  to  the  spheroid’s  symmetry.  By  varying  the  orientation  9  we  are  able  to  trace  the 
change  in  dipole  magnitude  and  angle  with  respect  to  the  applied  field.  Figure  4  shows 
these  curves  as  polar  plots  for  five  different  ordnance  items;  a  60-mm  mortar,  a  76-mm 
projectile,  an  81-mm  mortar,  a  105-mm  projectile  and  a  155-mm  projectile.  Fitting 
a  dipole  to  observed  data  will  produce  a  single  point  in  this  polar  plot  that  needs  to 
be  matched  to  an  ordnance  item.  The  curves  for  the  76-mm  projectile  and  the  81- 
mm  mortar  are  very  similar  in  certain  parts  of  the  plot,  indicating  that  discriminating 
between  these  items  may  be  very  difficult.  Further,  certain  orientations  lead  to  identical 
dipole  moments  of  both  these  items  to  the  60-mm  mortar  and  the  105-mm  projectile. 

The  analysis  in  the  last  paragraph  indicates  that  there  can  be  ambiguity  in  identi¬ 
fication  using  magnetometry.  This  fact  is  rigorously  established  in  Appendix  C  where 
it  is  shown  that  for  a  given  spheroid  at  a  particular  orientation,  there  are  an  infinite 
number  of  other  spheroids  that  could  have  produced  the  same  dipole  moment.  For  ex¬ 
ample,  Figure  5  and  Table  1  show  the  family  of  spheroids  that  can  produce  the  same 
dipole  moment  as  a  105-mm  projectile  orientated  at  45°  to  the  Earth’s  field.  The  90-mm 
projectile  lies  very  close  to  this  curve  emphasizing  the  difficulty  in  separating  these  two 
items.  The  folding  over  of  this  curve  to  encompass  spheroids  of  larger  dimension  is  of 
more  concern  to  classification;  the  155-mm  lies  quite  close  to  this  branch  of  the  curve. 
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There  is  also  a  family  of  oblate  spheroids  that  could  produce  the  same  dipole  moment; 
as  could  objects  of  different  shapes.  This  inherent  ambiguity  has  a  significant  impact 
on  the  methods  we  develop  for  ordnance  classification  in  the  remainder  of  this  report. 


Figure  5:  Spheroid  dimensions  that  can  produce  the  same  dipole  as  a  105-mm  shell  at 
45°  inclination  to  the  Earth’s  field.  The  dimensions  and  orientation  angles  corresponding 
to  the  labelled  points  are  shown  in  Table  1. 


Diameter  (mm) 

Aspect  ratio 

Angle  (degrees) 

82 

5 

40.5 

138 

5 

84.0 

138 

2.82 

57.2 

D 

327 

0.06 

53.4 

E 

327 

0.31 

33.4 

Table  1:  Spheroid  dimensions,  and  their  angles  relative  to  the  Earth’s  field,  that  produce 
the  same  dipole  moment  as  a  105-mm  projectile  at  45°  inclination. 
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3  Inversion  method  for  parameter  estimation 


The  previous  section  demonstrated  that,  in  many  situations,  it  would  be  difficult  to 
recover  more  than  the  dipole  component  of  a  spheroid’s  magnetic  anomaly.  Therefore, 
rather  than  attempting  a  direct  inversion  for  the  spheroid  location  and  dimension,  we 
apply  a  two-step  procedure: 

1.  Invert  for  the  best-fitting  dipole  moment  and  location; 

2.  Use  the  recovered  dipole  moment  for  discrimination  and  ordnance  identification. 

Due  to  possible  difficulties  in  properly  removing  the  regional  and  local  geologically 
derived  fields  from  a  magnetic  survey,  we  also  allow  for  a  dc-shift  in  our  dipole  model; 
this  means  that  seven  parameters  are  required  to  characterize  each  anomaly. 


Inversion  goal:  Find  the  location  x  =  ( x,y,z ),  dipole  moment  m  =  (rnx,my, mz) 
and  dc-shift  d,  that  produces  the  best  fit  to  a  set  of  N  total-field  observations, 
{dobs(xj),i  —  1, . . .  ,N}. _ 


We  follow  the  usual  convention  in  inverse  problems  and  use  m  to  specify  the  seven 
model  parameters  (it  should  be  clear  from  the  context  whether  m  refers  to  the  model 
parameters  or  the  dipole  moment).  We  formulate  the  inverse  problem  as  a  weighted, 
least-squares  optimization  procedure  where  the  goal  is  to  solve 

min  <£(m)  (12) 

me5?"  v  ’ 

The  objective  function  is 

4>(m)  =  ^(F(m)  -  do6s)TW(F(m)  -  d0*5)  (13) 

where  F(m)  is  the  forward  model  that  produces  the  predicted  data,  do6s  is  the  observed 
data  and  W  is  the  data-weighting  matrix  (ideally  the  data  weighting  matrix  is  the  inverse 
covariance  matrix  of  the  errors).  We  assume  independent,  gaussian  errors,  so  that  the 
data-weighting  matrix  is  diagonal  with  entries  Wu  =  1/of.  The  objective  function  may 
then  be  rewritten  as 


4>{m) 


1  N 


Fi{ m)  -  d 


obs 


T  2 


1  N 

^I>i(m)2 


(14) 


i=l 


where  rj(m)  is  the  i-th  normalized  residual,  and  we  are  assuming  that  there  are  N 
observations. 
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Inversion  method 

We  use  the  interior-reflective  Newton  method  (Branch  et  al.,  1999)  to  find  the  best 
fitting  model  parameters.  The  algorithm  is  a  trust  region  method  which  means  that  the 
vector  Sfc  that  updates  the  current  model,  i.e., 


=  mfe  +  sk  (15) 

is  chosen  from  a  quadratic  model  of  the  objective  function  within  a  region  where  this 
model  can  be  trusted  (hence  the  term  trust  region  method).  A  quadratic  approximation 
to  0(m)  is 

<Kmfc  +  sk)  «  <f>(mk)  +  V<t>(mk)sk  +  ^sk  V2(j)( mk)sk  (16) 

For  a  least  squares  problem  we  have 


V0(m)  =  JT(m)r(m)  (17) 

where  J  is  the  Jacobian  matrix  which  expresses  the  gradient  of  the  residuals  with  respect 
to  each  model  parameter, 


•Mm)  = 


dn(  m) 
dnij 


and  the  Hessian  matrix,  H,  is  given  by 


(18) 


N 

H(  m)  =  V2<£(m)  =  JT(m)  J(m)  +  ^  rj(m)V2rj(m)  (19) 

t=i 

The  Hessian  can  be  difficult  to  calculate  because  it  requires  second-order  derivatives. 
In  least-squares  problems  a  Gauss-Newton  approximation  is  often  used,  whereby  the 
second  term  in  the  above  equation  is  ignored.  This  is  a  good  approximation  when  the 
residuals  are  small  and  means  that  the  Hessian  can  be  obtained  essentially  for  free  once 
the  Jacobian  has  been  calculated. 

The  trust  region  method  involves  finding  the  s  (we  drop  the  subscripts  denoting 
iteration  number  and  the  dependence  on  m  for  clarity),  that  satisfies 

min  |^sTJTJs  +  sTJTr|  such  that  ||Ds||  <  A  (20) 

where  D  is  a  diagonal  scaling  matrix  and  A  is  a  positive  scalar  that  controls  the  trust- 
region  size.  Without  the  trust  region  constraint  the  solution  of  the  above  equation  can 
be  found  from  the  normal  equations 


JTJs  =  -JTr  (21) 

The  usual  trust-region  approach  is  the  Levenburg-Marquart  algorithm  (Marquardt, 
1963)  where  we  solve 


(. JTJ  +  AJ)s  =  -  JTr  (22) 

The  positive  parameter  A  ensures  that  the  solution  meets  the  constraint  \\Ds\\  <  A  and 
also  improves  the  condition  number  of  the  normal  equations. 
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Instead  of  minimizing  Equation  (20)  over  the  full  space,  the  algorithm  we  use  min¬ 
imizes  over  a  two-dimensional  subspace,  span(si,S2).  The  first  basis  of  the  subspace  is 
chosen  as  the  steepest  descent  direction 

si  =  -JT  r  (23) 

and  the  second  is  obtained  by  applying  preconditioned  conjugate  gradients  (PCG)  to 
Equation  (21).  The  PCG  iterations  are  terminated  when  a  negative  curvature  direction 
is  discovered,  i.e., 


sj  JTJs2  <  0  (24) 

or  until  the  residual  is  sufficiently  small;  this  would  be  an  approximate  Gauss-Newton 
step.  The  philosophy  behind  the  method  is  to  force  global  convergence  (by  the  steep¬ 
est  descent  or  negative  curvature)  while  still  achieving  rapid  local  convergence  (by  the 
Gauss-Newton  step). 

Bounds  on  model  parameters 

The  Interior-reflective  Newton  Method  allows  bounds  to  be  placed  on  the  model  pa¬ 
rameters  by  the  user  if  they  are  required.  We  typically  impose  bounds  on  the  spatial 
location  of  the  dipole,  e.g. 

-2  <  x  <  2,  -2  <  y  <  2  and  -  4  <  z  <  0  (25) 

but  do  not  bound  the  dipole  moment  or  dc-shift.  With  lower  and  upper  bounds  of  [lt,  iq], 
the  minimization  problem  of  Equation  (12)  is  changed  to 

min  (0(m)  :  L  <  m;  <  (26) 

m eft"  i  -  i  -  u  \  i 

This  constrained  minimization  problem  can  be  solved  using  a  modification  of  the  for¬ 
mulation  presented  above  (Branch  et  al.,  1999). 

Starting  model 

Due  to  the  possible  presence  of  local  minima  it  is  important  to  obtain  a  good  first  guess 
of  the  model  parameters.  To  achieve  this  goal,  we  use  the  approach  of  McFee  and  Das 
(1986),  whereby  the  dipole  parameters  are  analytically  obtained  from  the  z-component 
of  the  vector  field.  Typically,  the  data  available  comprise  the  total-field  and  we  first  need 
to  estimate  the  z-component.  Most  anomalies  from  ordnance  have  peak  amplitudes  of 
less  than  1,000  nT,  whereas  the  Earth’s  field  is  on  the  order  of  50,000  nT.  The  anomalous 
total  field  A Bt  is  then  approximately, 


A Bt  ~  axABx  +  ayABy  +  azABz  (27) 

where  a  =  (ax,ay,az)  is  a  unit  vector  in  the  direction  of  the  Earth’s  field  and  AB  = 
(ABX,  ABy,  ABZ)  is  the  vector  magnetic  anomaly.  Almost  all  of  the  continental  USA 
and  Canada  has  inclinations  exceeding  60°  so  that  az  >  0.85  and  to  a  first  approximation 

A  Bz~^  (28) 

In  regions  where  the  z-component  of  the  field  is  not  so  dominant,  A Bz  can  be  obtained 
from  A Bt  using  well  known  Fourier  transformation  techniques  (Blakely,  1996). 
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The  McFee  and  Das  (1986)  technique  involves  first  locating  the  maximum  and  min¬ 
imum  values  of  the  field.  Then 

•  The  fine  joining  these  two  points  gives  the  azimuth  of  the  dipole; 

•  The  distance  between  the  two  points  is  related  to  the  depth  of  the  dipole; 

•  The  locations  of  the  extremum  and  their  relative  amplitude  can  be  used  to  infer 
the  horizontal  location  (assumed  to  lie  on  the  fine  joining  the  points); 

•  The  relative  amplitudes  of  the  extremum  are  related  to  the  dipole  dip  angle; 

•  The  maximum  anomaly  amplitude  is  proportional  to  the  dipole  magnitude;  and 

•  Comparing  the  fit  of  the  dipole  to  the  observed  data,  estimates  the  dc-shift. 

For  a  nearly  vertical  dipole  there  is  only  a  well  defined  maximum  (or  minimum). 
In  that  case,  the  maximum  gives  the  horizontal  location,  the  half-width  relates  to  the 
depth  and  the  maximum  amplitude  scales  with  the  magnitude  of  the  vertical  dipole. 

Data  weighting  matrix 

For  the  least-squares  data  fitting,  we  weight  each  data-point  d°bs,  with  an  estimate  of 
the  standard  deviation  at  that  point,  cq,  through  the  data-weighting  matrix, 

W-  =  ±  = _ I _ 

of  (e  +  j\dfs\)2 

The  parameter  7  is  a  percentage  error  that  down-weights  the  influence  of  the  large  data 
values;  typically  we  chose  1%  <  7  <  5%.  The  parameter  e  ensures  that  the  denominator 
in  the  above  equation  does  not  get  too  close  to  zero  when  approaches  zero;  typically 
we  choose  e  InT. 

Model  parameter  scaling 

If  different  model  parameters  have  significantly  different  scale  lengths  and  typical  values 
(e.g.  position  ~  1  m;  dipole  moment  ~  0.01  Am2),  then  small  favorable  changes  in  some 
model  parameters  may  be  hidden  by  the  large  changes  required  in  others.  The  ideal  is  to 
have  a  non-dimensionalized  model  space,  where  unit  changes  in  each  model  parameter 
have  about  the  same  impact  on  the  objective  function.  We  therefore  scale  each  model 
parameter  by  dividing  by  the  following  values, 

xij/p  =  (1, 1,  l)m,  m typ  =  (0.1,0.1,0.1)Am2  and  dtyp  =  5  nT  (30) 
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4  Application  of  the  inversion  method 


Currently,  the  inversion  algorithms  are  implemented  in  Matlab  and  comprise  part  of 
a  GUI  driven  toolbox  for  geophysical  processing  developed  at  The  University  of  British 
Columbia.  Anomalies  are  first  identified  in  a  data-set  using  a  combination  of  automatic 
and  manual  interpretation.  Then,  the  inversion  algorithm  is  applied  to  a  user  specified 
area  about  each  anomaly  (typically  around  2  x  2  m2  to  3  x  3  m2).  The  algorithm  is 
terminated  when  the  objective  function  decreases  by  less  than  a  specified  tolerance, 
or  when  the  gradient  falls  below  a  specified  value.  Both  these  parameters  are  user 
controlled. 

Summary  reports  can  be  generated,  and  the  results  visually  inspected  as  demon¬ 
strated  in  Figure  6.  We  will  illustrate  the  performance  of  the  inversion  algorithm  on  a 
large  magnetics  survey  collected  as  part  of  an  clearance  project  in  the  Helana  Valley, 
Montana. 


Raw  data  575 


Model  tit  575 


Residual  575 


X  =  33.43  (±0.010)  m 
Y=  10.71  (±0.010)  m 
Z  =  -0.25  (±0.008)  m 
Moment  =  0.0514  (±0.0017)  Am2 
Azimuth  =  -98.6°  (±6.1°) 

Dip  =  73.2°  (±1.7°) 

CorrCoeff  =  0.970 
Angle  =  152.8°  (±1.7°) 

Successful  Fit 


Figure  6:  Example  result  of  dipole  moment  inversion;  the  numbers  in  brackets  are  the 
estimated  standard  deviations  of  the  recovered  parameters.  Azimuth  and  Dip  refer  to 
the  orientation  of  the  dipole  in  geographic  coordinates,  while  Angle  is  measured  between 
the  dipole  moment  and  the  Earth’s  field. 
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Figure  7:  Some  of  the  81-mm  mortars  recovered  at  Guthrie  Rd. 


The  Guthrie  Road,  Montana  dataset 

Once  open  fields,  portions  of  the  Helena  Valley  were  used  in  the  1950s  for  military 
training  by  Montana  Army  National  Guard  (MTARNG).  The  National  Guard  Bureau 
funded  the  remediation  project  which  was  divided  into  two  phases  focusing  on  different 
impact  areas:  Diamond  Springs  and  Guthrie  Road.  We  consider  the  data  collected  in 
the  Guthrie  Road  area  in  the  summer  of  1998.  G-tek  Australia,  conducted  a  magne¬ 
tometer  survey  of  the  area  using  an  all-terrain  vehicle  towing  an  array  of  eight  cesium 
vapor  magnetometers  (Clark  et  al.,  1999).  The  system  was  equipped  with  a  real-time 
global  positioning  system  that  allowed  anomalies  to  be  positioned  to  within  about  20 
centimeters. 

G-tek  found  840  anomalies  that  were  tagged  as  potential  UXO,  and  these  were  ex¬ 
cavated  in  the  summer  of  1999.  We  analyze  822  of  the  anomalies  (we  couldn’t  locate 
the  data  for  the  remaining  18  anomalies).  Table  2  partitions  the  dig-sheet  identifica¬ 
tions  into  six  categories;  76-mm  projectiles,  81-mm  mortars,  large  pieces  of  ordnance, 
shrapnel,  metallic  debris  and  geology  (including  anomalies  where  no  metallic  items  were 
found). 

A  selection  of  the  81-mm  mortars  recovered  at  Guthrie  Rd  are  shown  in  Figure  7. 
Most  are  intact  and  have  not  been  substantially  deformed  by  the  force  of  impact.  How¬ 
ever,  there  are  a  number  of  deformed  but  intact  items.  These  tend  to  have  been  practice 
rounds  that  were  deformed  when  the  black  spotter  charge  detonated  on  impact.  There 
are  also  a  number  of  items  that  have  been  split  into  pieces,  these  were  not  classified  as 
UXO. 

We  fitted  dipole  moments  using  a  percentage  error  of  7  =  5%  and  a  base  error  of 
e  =  1  nT.  The  entire  run  of  inversions  took  around  |  hour  on  a  fairly  standard  desktop 
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computer  (Pentium  Pro  Processor  with  800  MHz  processor  and  1  Gb  RAM).  We  use 
the  correlation  coefficient  c  between  fitted  and  observed  data  to  determine  the  goodness 
of  fit.  We  found  that  85,  or  around  10%  of  anomalies  had  fits  with  c  <  0.7  and  162,  or 
around  20%  had  c  <  0.8.  A  break-down  according  to  the  dig-sheets  of  the  failed  fits  is 
given  in  Table  2.  Most  of  the  failed  fits  are  from  items  identified  as  shrapnel,  geology 
or  metallic  debris,  with  only  one  ordnance  item  having  a  failed  fit.  We  did  not  visually 
inspect  every  failure,  but  those  we  did  inspect  generally  had  some  type  of  data  issue. 
These  included  overlapping  acquisition  lines  that  we  couldn’t  remove  without  undue 
effort;  levelling  problems  between  adjacent  lines;  and  overlapping  anomalies.  In  no  case 
did  the  inversion  algorithm  fail  where  the  data  quality  was  good.  In  the  next  sections 
we  will  describe  methods  for  discrimination  and  identification  and  apply  them  to  this 
data-set. 


76-mm 

81-mm 

Large  pieces 

Shrapnel 

Metal 

Geology 

Total 

All  items 

34 

49 

19 

217 

377 

126 

822 

Failed  fits 

0 

1 

1 

13 

41 

29 

85 

Table  2:  Number  of  items  found  during  excavation,  and  the  number  of  fits  with  c  <  0.7 
(failed  fits) 
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5  UXO  Discrimination 


The  last  section  described  how  dipole  moments  and  position  are  fit  to  each  anomaly. 
If  the  data-fit  is  unsatisfactory  the  anomaly  is  marked  as  such  and  in  a  production  clean¬ 
up  would  be  excavated  as  potential  UXO.  Prom  the  viewpoint  of  data  interpreters,  we 
are  unable  to  provide  any  reliable  information  on  the  target  characteristics.  On  the 
other  hand,  if  the  fit  is  satisfactory  we  attempt  to  use  the  recovered  dipole  moment 
to  discriminate  UXO’s  from  non-UXO’s.  For  those  items  classified  as  potential  UXO 
the  additional  step  of  ordnance  identification  is  attempted  as  described  in  the  following 
section. 

Our  method  for  ordnance  discrimination  is  based  on  that  described  earlier  in  this 
report,  and  is  due  to  Nelson  et  al.  (1998).  Briefly,  the  method  relies  on  the  following 
two  factors; 

1.  Shock  demagnetization  causes  the  remnant  magnetism  of  an  ordnance  to  be  vir¬ 
tually  eliminated; 

2.  The  induced  magnetization  in  a  spheroid  with  e  <  7  is  constrained  to  he  within 
about  60°  of  the  Earth’s  field  (Figure  3). 

With  a  recovered  dipole  moment  of  m,  our  discrimination  routine  is  to  calculate  the 
angle  between  it  and  the  Earth’s  field,  Bc, 

* = arccos(i5rEii)  <31) 

We  can  never  be  certain  that  all  remnant  magnetization  has  been  erased  by  the  shock 
of  impact.  Additionally,  non-ordnance  items,  whether  remnantly  magnetized  or  not, 
may  have  dipole  angles  close  to  the  Earth’s  field.  Consequently,  it  is  not  possible  to 
make  a  definitive  ordnance/non-ordnance  distinction  for  all  anomalies.  Our  preferred 
discrimination  method  is  to  rank  items  according  to  the  likelihood  they  are  due  to 
ordnance.  We  do  this  on  the  basis  of  the  angle  ip]  the  items  with  smaller  ip  are  ranked 
higher  up  the  list.  In  a  clean-up  operation  one  then  excavates  according  to  the  list.  If 
resources  are  available  to  dig  50%  of  the  anomalies  then  the  first  50%  of  items  in  the 
fist  are  excavated.  As  a  quality  control  measure  a  certain  number  of  the  lower  ranked 
items  should  also  be  excavated. 

The  Guthrie  Road,  Montana  dataset 

In  Figure  (8a)  we  plot  the  dipole  magnitude  and  angle  relative  to  the  Earth’s  field  of  all 
the  non-ordnance  anomalies.  Most  of  the  dipoles  he  within  ±90°  of  the  Earth’s  field  but 
a  large  number  have  angles  exceeding  60°.  On  the  other  hand,  the  dipole  angles  for  the 
ordnance  and  large  ordnance  fragments  almost  all  he  within  ±60°  of  the  Earth’s  field 
Figure  (8b).  These  results  indicate  that  significant  shock  demagnetization  has  occurred 
to  the  ordnance.  There  is  only  one  significant  outlier,  with  an  angle  of  around  70°. 

In  Figure  (9a)  we  give  a  normalized  histogram  of  the  angles  for  the  different  categories 
of  items.  There  is  a  clear  concentration  of  ordnance  with  dipole  angles  close  to  the 
Earth’s  field  (we  include  the  significant  sized  ordnance  fragments  in  this  category). 
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Figure  8:  (a)  Polar  plot  of  dipole  angle  and  magnitude  (Am2)  of  non-ordnance  items; 
(b)  Same  for  ordnance  (cross)  and  significant  ordnance  fragments  (dots). 


20  40  60  80  100  120  0 

Dipole  angle  relative  to  Earth'  field 


■■  Ordnance 

Shrapnel 

■i  Metallic  debris 
I  I  Geology 


0.1  0.2  0.3  0.4 

Dipole  magnitude  (Am2) 


Figure  9:  Breakdown  of  dig-sheets  identifications  and  their  (a)  dipole  angle  relative  to 
the  Earth’s  magnetic  field;  and  (b)  dipole  magntitude. 


Shrapnel  tends  to  have  relatively  low  angles,  but  there  are  a  significant  number  of  items 
with  larger  angles.  There  are  two  simple  mechanisms  that  could  explain  this  distribution: 

1.  Shock  demagnetization  sometimes  occurs  with  shrapnel  but  not  always; 

2.  Shrapnel  acquires  a  remnant  magnetization  after  impact,  in  the  presence  of  the 
Earth’s  field. 

Metallic  debris  and  geological  anomalies  can  have  significant  remnant  magnetism, 
with  a  definite  peak  in  the  distribution  between  60°  and  100°.  We  can  think  of  no 
plausible  explanation  for  this  peak  in  the  distribution  of  dipole  angles. 
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Angle 

Dig 

Leave 

Ordnance  found 

Ordnance  missed 

False  alarms 

449 

373 

78 

5 

371 

485 

337 

81 

2 

404 

65 

522 

300 

82 

1 

440 

70 

553 

269 

82 

1 

471 

75 

605 

217 

83 

0 

522 

Table  3:  Guthrie  Road  discrimination  results  for  intact  ordnance  items.  Discrimination 
used  the  dipole  angle  only. 


Angle 

Dig 

Leave 

Ordnance  found 

Ordnance  missed 

False  alarms 

449 

373 

91 

11 

358 

485 

337 

99 

3 

386 

65 

522 

101 

1 

70 

553 

269 

101 

1 

452 

75 

605 

217 

102 

0 

503 

Table  4:  Guthrie  Road  discrimination  results,  including  large  pieces  of  ordnance  items. 
Discrimination  used  the  dipole  angle  only. 


The  results  of  discrimination  using  different  cut-off  angles,  for  intact  ordnance  items 
are  summarized  in  Table  3,  and  for  ordnance  items  and  large  fragments  in  Table  4.  The 
tables  give  the  number  of  items  to  dig  (including  the  85  anomalies  with  poor  fits),  the 
number  left  in  the  ground,  the  ordnance  found,  the  ordnance  left  and  the  number  of 
false  alarms.  Figure  (10)  plots  the  proportion  of  items  recovered  versus  the  number  of 
holes  that  need  to  be  excavated.  Once  the  85  anomalies  with  poor  fits  have  been  dug, 
there  is  a  high  yield  of  ordnance  as  the  number  of  holes  increases.  However,  once  about 
80%  of  items  have  been  recovered  there  is  a  noticeable  slow-down  in  the  ordnance  yield; 
a  lot  of  holes  need  to  be  dug  to  recover  the  last  few  ordnance  items.  A  total  of  560 
(68%)  of  holes  are  required  to  recover  all  ordnance  items;  this  number  also  recovers  all 
the  large  ordnance  pieces. 

Many  of  the  dipole  fits  for  shrapnel  and  metallic  debris  have  moments  with  low 
magnitude  (Figure  9b).  None  of  the  intact  ordnance  items  has  a  dipole  magnitude  less 
than  0.05  Am2.  However,  seven  of  the  large  ordnance  fragments  have  dipole  moments 
less  than  this  value.  The  results  of  rejecting  items  as  non-ordnance  based  on  the  dipole 
magnitude  (as  well  as  a  cutoff  angle  of  75°)  are  summarized  in  Tables  5  (for  ordnance 
only)  and  6  (for  large  fragments  as  well).  With  a  dipole  cutoff  of  0.05  Am2  no  intact 
ordnance  items  are  missed  and  we  can  leave  379  items  in  the  ground  (46%  saving).  If 
we  make  the  reasonable  assumption  (for  this  area)  that  the  smallest  item  likely  to  be 
found  is  a  60-mm  mortar,  the  0.05  Am2  cutoff  is  quite  sensible;  the  minimum  dipole 
moment  expected  from  a  60-mm  mortar  is  around  0.055  Am2.  Ranking  items  by  their 
magnitude  results  in  a  lower  initial  yield  of  ordnance  (Figure  10)  than  ranking  them  by 
angle.  However,  the  magnitude  ranking  does  not  suffer  as  significant  a  slow-down  in 
yield,  with  a  crossover  occurring  once  about  80%  of  items  have  been  recovered.  A  total 
of  443  (54%)  holes  need  to  be  excavated  to  recover  all  the  ordnance. 

One  additional  discrimination  method  is  examined.  It  uses  the  idea  of  remnant 
magnetism  and  requires  information  on  the  ordnance  items  likely  to  be  found  in  the 
area.  For  the  Guthrie  Rd  data  we  assume  that  60  and  81-mm  mortars  and  76,  90,  105 
and  155-mm  projectiles  may  be  present  in  the  area  and  we  use  the  ordnance  dimensions 
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Figure  10:  Yield  of  ordnance  with  number  of  holes  dug  for  the  three  discrimination 
methods.  For  all  methods  it  is  assumed  the  85  anomalies  with  poor  fits  are  excavated 
first. 


Mag.  (A/m) 

Dig 

Leave 

Ordnance  found 

Ordnance  missed 

False  alarms  | 

'  ,  ■ 

383 

.  . 

81 

2 

443 

83 

0 

Table  5:  Guthrie  Road  discrimination  results  for  intact  ordnance  items.  Discrimination 
used  the  dipole  angle  (75°)  and  magnitude. 


fisted  in  Table  7.  The  method  is  described  in  more  detail  in  the  next  section.  Briefly, 
we  determine  the  minimum  percentage  of  remnant  magnetization  required  to  make  the 
recovered  dipole  moment  match  of  the  ordnance  items.  We  then  rank  the  items  according 
to  the  percentage  of  remnant  magnetism,  with  lower  percentages  higher  up  the  fist.  The 
method  does  not  apply  to  large  ordnance  fragments  as  these  are  not  considered  in  our 
modelling. 

The  results  of  discriminating  by  remnant  magnetism  are  summarized  in  Table  8. 
Note,  that  we  also  incorporate  the  angle  constraint  (<  75°)  and  the  magnitude  con¬ 
straint  (>  0.05  Am2).  Allowing  up  to  half  of  the  dipole  anomaly  to  be  due  to  remnant 
magnetization,  we  can  identify  all  ordnance  items  and  leave  420  anomalies  in  the  ground 
(51%  saving).  The  ordnance  versus  number  of  holes  curve  for  remnant  magnetization 
(Figure  10)  reveals  a  high  yield  of  ordnance  with  increasing  number  of  excavations. 
There  is  only  a  very  slight  slowdown  in  yield  once  about  95%  of  ordnance  items  have 
been  found.  Once  402  (49%)  of  anomalies  are  excavated  all  ordnance  items  are  recovered. 
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Mag.  (Am2) 

Dig 

Leave 

Ordnance  found 

Ordnance  missed 

False  alarms 

0.06 

383 

439 

12 

293 

0.05 

443 

379 

95 

7 

348 

0.04 

493 

329 

98 

4 

395 

0.03 

558 

264 

100 

2 

458 

0.02 

597 

225 

102 

0 

495 

Table  6:  Guthrie  Road  discrimination  results,  including  large  pieces  of  ordnance  items. 
Discrimination  used  the  dipole  angle  (75°)  and  magnitude. 


Item 

Diameter  (mm) 

Aspect  ratio 

60-mm  mortar 

60 

3.78 

76-mm  projectile 

76 

3.57 

81-mm  mortar 

81 

3.54 

90-mm  projectile 

90 

5.0 

105-mm  projectile 

105 

3.71 

155-mm  projectile 

155 

4.0 

Table  7:  Ordnance  sizes  used  for  the  Guthrie  Road  identification  results. 


Remnant  (%) 

Dig 

Leave 

Ordnance  found 

Ordnance  missed 

False  alarms 

10 

228 

594 

46 

37 

182 

20 

306 

516 

71 

12 

235 

30 

350 

472 

79 

4 

271 

40 

384 

438 

82 

1 

302 

50 

415 

407 

83 

0 

332 

Table  8:  Guthrie  Road  discrimination  results  for  intact  ordnance  items.  Discrimination 
used  the  dipole  angle  (75°),  magnitude  (<  0.05  Am2)  and  remnant  magnetism. 
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6  Ordnance  identfication 


Identification  using  the  recovered  dipole  moment 

For  those  anomalies  classified  as  potential  ordnance,  the  next  step  is  to  try  to  identify 
the  ordnance  type.  We  established  in  Appendix  B  and  an  earlier  part  of  the  report  that 
a  given  dipole  anomaly  could  have  been  produced  by  an  infinite  number  of  spheroids 
(Figure  5).  However,  for  any  given  site,  there  are  only  a  finite  number  of  ordnance 
types  that  are  likely  to  be  found.  This  is  the  key  observation  that  allows  us  to  identify 
ordnance  type  in  the  presence  of  ambiguity. 

The  dipole  induced  in  a  spheroid  depends  on  the  angle  between  the  spheroids  semi¬ 
major  axis  and  the  Earth’s  field.  There  is  no  azimuthal  dependence  due  to  symmetry. 
As  the  ordnance  is  rotated  about  the  Earth’s  field  the  induced  dipole  moment  will  trace 
out  a  curve  such  as  that  shown  in  Figure  (4).  This  represents  all  dipole  moments  that 
can  be  produced  by  the  given  ordnance  item  in  the  absence  of  remnant  magnetism. 
This  mapping  procedure  can  be  repeated  for  all  ordnance  items  suspected  of  occurring 
in  the  area.  With  N0  ordnance  items,  a  series  of  curves,  k  =  1, . .. ,  N0},  will 

be  generated.  For  a  particular  recovered  dipole  moment  m,  we  calculate  the  minimum 
distance  between  it  and  each  of  the  dipole  curves 

Amfc=  min  ||m-mfc(0)||  (32) 

O<0<27r 

We  assume  that  any  discrepancy  is  due  to  remnant  magnetism  which,  as  a  percentage 
of  the  observed  moment,  will  be 


Ik  =  100 


A  mk 

IMI 


(33) 


This  allows  us  to  rank  the  ordnance  items  according  to  the  likelihood  they  produced 
the  observed  anomaly.  The  most  likely  item  is  that  with  the  lowest  percentage  remnant 
magnetism. 

The  success  of  the  procedure  will  depend  on; 


1.  The  accuracy  of  the  recovered  dipole  moment  (see  the  error  analysis  section); 

2.  The  orientation  of  the  ordnance,  because  some  orientations  produce  dipoles  that 
are  close  to  dipoles  potentially  produced  by  other  ordnance  items; 

3.  How  closely  the  spheroid  approximation  matches  that  for  real  ordnance;  and 

4.  How  much  remnant  magnetism  the  ordnance  maintains  after  the  shock  of  impact. 


The  Guthrie  Road,  Montana  dataset 

The  dipole  moments  for  those  anomalies  due  to  76-mm  projectiles  are  shown  in  Fig¬ 
ure  (11a)  and  for  81-mm  mortars  in  Figure  (lib).  Misidentifications  for  both  types  of 
ordnance  appear  likely  due  to  the  overlapping  curves  for  the  difference  ordnance  types. 
In  particular  the  moments  for  almost  all  the  76-mm  projectiles  are  in  a  region  where  the 


22 


76  81  90  105 

Ordnance  diameter  (mm) 


76  81  90  105 

Ordnance  diameter  (mm) 


Figure  12:  Identification  results  for  the  Guthrie  Road  data  using  both  the  dipole  and 
dipole  plus  octupole  methods  for  (a)  76-mm  projectiles;  and  (b)  81-mm  mortars. 


Dipole  identification  results  for  76-mm  projectiles  are  given  in  (Figure  12a)  and  for 
81-mm  mortars  in  Figure  (12b).  In  19  out  of  34  cases  (56%)  the  76-mm  projectile  is 
identified  correctly.  The  most  common  case  of  misidentilication  is  as  a  60-mm  mor¬ 
tar;  from  Figure  (4)  we  see  that  this  is  due  to  overlapping  moments  between  the  60 
and  76-mm  ordnances.  For  81-mm  mortars  the  identification  routine  is  correct  in  25 
of  48  cases  (52%),  so  that  overall  we  achieve  correct  identifications  54%  of  the  time. 
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Figure  13:  One  of  the  practice  81-mm  mortars  deformed  due  to  the  detonation  of  the 
spotting  charge. 

The  most  common  misidentification  for  the  81-mm  mortar  is  as  a  76-mm  projectile  (13 
times).  Some  of  the  errors  in  misidentification  could  relate  to  the  deformation  of  cer¬ 
tain  of  the  81-mm  practice  mortars  that  occurred  when  the  spotter  charge  detonated 
(Figure  13).  A  spheroid  is  unlikely  to  give  a  very  good  approximation  to  this  item.  Un¬ 
fortunately,  while  recovered  ordnance  items  were  retained  after  excavation  no  records 
of  the  anomaly  number  were  kept.  Thus,  we  cannot  determine  which  recovered  dipole 
moments  correspond  to  deformed  ordnance. 

Assuming  we  are  using  the  remnant  magnetization  method  for  discrimination  and  use 
a  cutoff  of  50%,  there  will  be  235  false  alarms;  items  incorrectly  identified  as  ordnance. 
Applying  the  identification  method  to  these  false  alarms  (Figure  14a)  is  most  likely  to 
cause  a  misidentification  as  a  60-mm  mortar  (the  moments  are  generally  quite  low). 
Metallic  debris  is  more  likely  to  give  a  false  identification  as  a  large  ordnance  item  than 
any  of  the  other  categories  of  items  (shrapnel,  large  ordnance  pieces  and  geology). 

Seeded  test  site  at  The  Former  Fort  Ord,  California 

The  former  Fort  Ord  is  located  near  Monterey  Bay  in  northwestern  Monterey  County, 
California.  Since  1917,  portions  of  Fort  Ord  were  used  by  Army  units  for  maneuvers, 
target  ranges,  and  other  purposes.  Investigation  and  ordnance  removal  actions  have 
been  conducted  at  the  former  Fort  Ord  since  1994  when  the  base  closed.  In  November 
1998,  the  United  States  Department  of  the  Army  instigated  an  Ordnance  Detection  and 
Discrimination  Study  (ODDS)  at  Fort  Ord  (Asch  and  Staes,  2001). 

The  purpose  of  the  Seeded  Test  was  to  evaluate  the  detection  and  discrimination 
capabilities  of  selected  geophysical  instruments  when  ordnance  items  are  buried  or 
’’seeded”  below  the  ground  surface  in  a  local  site-specific  geologic  conditions.  The  seeded 
test  area  consists  of  five  separate,  adjacent  grids,  which  were  cleared  of  anomalies  to  pro- 
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Figure  14:  Identification  results  for  the  Guthrie  Road  data  for  (a)  dipole  method;  and 
(b)  dipole  plus  octupole  method. 


vide  a  controlled  testing  environment,.  Two  grids  are  ” known”  plots  where  the  items, 
depths  and  orientations  of  items  are  made  known  to  all  participants.  It  is  these  two 
known  sites  (K1  and  K2)  that  we  analyze  here. 

Seventeen  different  types  of  ordnance  items  were  buried  in  the  seeded  test  plots. 
Some  ordnance  scrap  and  other  items  were  also  seeded  to  allow  evaluation  of  discrimi¬ 
nation  capabilities.  Additional  areas  of  the  plots  were  excavated  and  back-filled  without 
any  buried  items. 

For  Kl,  5  of  the  7  ’’empty  hole”  items  (618,  619,  621  to  623)  did  not  produce 
acceptable  fits,  so  they  would  not  have  been  declared  as  potential  ordnance.  The  other 
two  ’’empty  holes”  had  poorly  fitting  dipoles  and  would  have  been  identified  as  scrap. 
A  fragment  (595)  and  a  signal  illumination  flare  (612)  also  produced  unacceptable  fits. 

Discrimination  and  identification  results  for  the  remaining  items  are  summarized  in 
Table  9.  Note,  that  the  column  marked  Interp  lists  the  first  choice  interpretation  for 
each  item  (e.g.  scrap,  22  mm  projectile,  60  mm  ordnance  etc).  The  column  marked 
Rank  shows  where  the  actual  ordnance  item  was  ranked.  For  instance,  consider  the 
2.36”  rockets.  If  the  rank  was  1  then  the  item  was  identified  correctly  as  a  2.36”  rocket 
(e.g.  anomaly  613).  If  the  rank  was  greater  than  1  (say  3),  then  the  2.36”  rocket  was 
ranked  as  the  third  most  likely  cause  of  the  anomaly  (e.g.  anomaly  602). 

Two  of  the  four  fragments  were  correctly  identified  as  non-ordnance  (using  a  discrim¬ 
ination  cutoff  of  60°).  Six  of  the  twelve  35-mm  projectiles  were  incorrectly  classified  as 
scrap;  this  may  be  due  to  remnant  magnetism.  Five  of  the  remaining  35-mm  projectiles 
had  the  correct  identification  ranked  as  one  of  the  top  three  choices.  One  of  the  2.36” 
rockets  was  classified  as  scrap,  the  other  three  had  a  2.36”  identification  ranked  in  the 
top  three  (one  was  correct).  The  six  signal  illumination  flares  could  not  be  identified 
correctly,  neither  could  any  of  the  grenades.  Overall,  the  discrimination/identification 
results  for  Kl  are  not  very  impressive.  We  believe  part  of  the  reason  could  be  due  to 
remnant  magnetism  (the  items  were  not  shock  demagnetized  prior  to  emplacement).  It 
may  also  be  a  reflection  of  the  generally  poor  performance  of  magnetics  for  small  calibre 
items  (Klaff  et  al.,  2002). 
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Number 

Item 

za  (m) 

zr  (m) 

c 

Angle 

Interp 

Rank 

620 

Empty  Hole 

0.00 

0.47 

0.79 

127.4 

Scrap 

617 

Empty  Hole 

0.00 

0.29 

0.75 

97.0 

Scrap 

599 

Fragment 

0.30 

0.92 

0.88 

86.1 

Scrap 

598 

Fragment 

0.30 

0.62 

0.90 

114.4 

Scrap 

597 

Fragment 

0.15 

0.28 

0.97 

48.8 

Signals 

596 

Fragment 

0.30 

0.41 

0.91 

55.2 

22  mm 

616 

35-mm 

0.61 

0.56 

0.93 

74.0 

Scrap 

608 

35-mm 

0.46 

0.23 

0.85 

59.0 

22  mm 

2 

607 

35-mm 

0.30 

0.23 

0.82 

150.1 

Scrap 

606 

35mm 

0.15 

0.25 

0.98 

33.6 

37-mm 

6 

605 

35-mm 

0.46 

0.95 

0.91 

52.4 

Signals 

3 

604 

35-mm 

0.61 

1.50 

0.81 

65.3 

Scrap 

603 

35-mm 

0.15 

0.37 

0.90 

55.2 

22  mm 

2 

600 

35-mm 

0.30 

0.43 

0.86 

68.8 

Scrap 

593 

35-mm 

0.61 

0.53 

0.92 

38.0 

22  mm 

3 

592 

35-mm 

0.30 

0.44 

0.97 

87.1 

Scrap 

591 

35-mm 

0.61 

0.34 

0.95 

105.6 

Scrap 

590 

35-mm 

0.46 

0.38 

0.96 

37.5 

60-mm  mortar 

2 

613 

2.36”  rocket 

0.61 

0.65 

0.96 

57.1 

Rocket  2  .36” 

1 

602 

2.36”  rocket 

0.46 

0.59 

0.98 

55.1 

Signals 

3 

594 

2.36”  rocket 

0.61 

0.89 

0.96 

50.8 

Rocket  3  .5” 

4 

586 

2.36”  rocket 

0.61 

2.79 

0.96 

118.2 

Scrap 

615 

Signals 

0.15 

0.26 

0.97 

39.7 

37-mm 

2 

614 

Signals 

0.30 

0.61 

0.82 

82.6 

Scrap 

611 

Signals 

0.15 

0.33 

0.89 

43.3 

22  mm 

3 

589 

Signals 

0.30 

1.08 

0.89 

20.1 

22  mm 

8 

588 

Signals 

0.30 

1.51 

0.95 

31.6 

81-mm  mortar 

6 

610 

Hand  grenade 

0.15 

0.08 

0.93 

49.9 

22  mm 

12 

609 

Hand  grenade 

0.30 

0.58 

0.77 

38.6 

22  mm 

12 

587 

Hand  grenade 

0.30 

0.88 

0.99 

121.5 

Scrap 

601 

Rifle  grenade 

0.46 

1.33 

0.78 

77.9 

Scrap 

Table  9:  Fort  Ord  identification  results  for  Kl,  where  za  is  the  actual  burial  depth,  zr 
is  the  recovered  depth,  c  is  the  correlation  coefficient  between  observed  and  predicted 
data  and  angle  is  the  angle  between  the  Earth’s  field  and  the  recovered  dipole  moment. 


For  the  K2  plot,  good  dipole  fits  could  not  be  obtained  for  four  items,  555,  556  (both, 
25-mm  projectiles  at  46  cm),  569  (fragment  at  30  cm)  and  574  (signals  at  30  cm).  The 
results  for  the  rest  of  the  items  are  summarized  in  Table  10.  There  were  seventeen  items 
correctly  identified  as  ordnance.  Six  of  these  were  identified  as  the  correct  ordnance 
type  and  five  others  had  the  correct  ordnance  listed  in  the  top  three  possibilities.  The 
other  six  identifications  were  a  long  way  off  and  we  believe  that  several  of  these  were 
remnantly  magnetized;  in  particular  all  the  stokes  mortars  and  the  hand  grenade. 

Identification  using  the  dipole  and  octupole 

Magnetic  anomalies  from  ordnance  are  usually  dominated  by  the  dipole  component; 
however,  with  high  quality  data  we  may  be  able  to  obtain  useable  information  from 
the  octupole  term.  Due  to  the  dipole  dominance  it  would  probably  not  be  a  good  idea 
to  attempt  a  direct  inversion  for  ordnance  dimensions  (local  minima,  non-uniqueness 
etc).  Rather,  a  better  method  is  to  obtain  an  optimum  fit  for  each  ordnance  item  in 
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Number 

Item 

Za  (m) 

zr  (m) 

c 

Angle 

Interp 

Rank 

565 

Fragment 

0.15 

0.14 

0.93 

92.1 

Scrap 

566 

Fragment 

0.30 

0.34 

0.98 

113.5 

Scrap 

567 

Fragment 

0.15 

0.14 

0.95 

52.3 

22  mm 

568 

Fragment 

0.30 

0.28 

0.89 

37.3 

37-mm 

578 

Fragment 

0.15 

0.22 

0.92 

32.5 

60-mm  mortar 

579 

Fragment 

0.30 

1.82 

0.86 

90.1 

Scrap 

581 

Fragment 

0.30 

0.57 

0.90 

36.4 

22  mm 

580 

25-mm 

0.61 

0.59 

0.82 

72.1 

Scrap 

553 

35-mm 

0.61 

0.67 

0.95 

29.1 

35-mm 

1 

559 

35-mm 

0.61 

0.77 

0.89 

33.2 

75-mm 

8 

575 

35-mm 

0.30 

0.26 

0.97 

154.0 

Scrap 

563 

37-mm 

0.46 

0.20 

0.81 

55.6 

22  mm 

3 

583 

37-mm 

0.61 

0.93 

0.85 

47.7 

35-mm 

5 

584 

37-mm 

0.61 

1.03 

0.89 

61.8 

Scrap 

554 

2.36”  rocket 

0.61 

0.89 

0.97 

21.6 

90-mm 

9 

558 

2.36”  rocket 

0.61 

0.81 

0.94 

49.5 

105-mm 

3 

564 

75-mm 

0.76 

0.67 

0.96 

12.6 

75-mm 

1 

572 

75-mm 

0.91 

1.85 

0.86 

85.6 

Scrap 

560 

Signals 

0.30 

0.50 

0.81 

124.4 

Scrap 

576 

Signals 

0.15 

0.16 

0.95 

99.8 

Scrap 

573 

3”  Stokes 

1.02 

0.63 

0.86 

60.0 

35-mm 

10 

577 

3”  Stokes 

1.02 

1.22 

0.89 

19.2 

155-mm 

2 

585 

3”  Stokes 

0.91 

0.66 

0.87 

16.2 

35-mm 

9 

561 

81 -mm 

0.91 

1.13 

0.96 

39.5 

81-mm  mortar 

1 

571 

81-mm 

0.91 

0.85 

0.94 

34.8 

90-mm 

2 

557 

90-mm 

0.76 

0.82 

0.97 

21.3 

90-mm 

1 

570 

90-mm 

0.76 

0.73 

0.97 

20.2 

90-mm 

1 

562 

105-mm 

1.22 

1.16 

0.94 

59.7 

Rocket  2.36” 

3 

582 

105-mm 

1.22 

1.50 

0.96 

17.8 

105-mm 

1 

552 

Hand  grenade 

0.30 

0.44 

0.92 

44.4 

Signals 

12 

Table  10:  Fort  Ord  identification  results  for  K2,  where  za  is  the  actual  burial  depth,  zr 
is  the  recovered  depth,  c  is  the  correlation  coefficient  between  observed  and  predicted 
data  and  angle  is  the  angle  between  the  Earth’s  field  and  the  recovered  dipole  moment 


our  library,  and  declare  the  ordnance  with  the  best  fit  as  the  most  likely  source  of  the 
anomaly.  For  each  item  in  the  library  this  procedure  involves  determination  of  six  model 
parameters  (x,  y,  z,  9,  <j>,  d),  where  (x,  y,  z)  is  location,  9  is  the  ordnance  dip  angle  relative 
to  horizontal,  <p  is  the  azimuth  angle  relative  to  geographic  North  and  d  is  the  dc-shift. 

We  again  use  a  least-squares  formulation  and  the  interior-reflective  Newton  method. 
The  data-weighting  and  parameter  bounds  are  identical  to  that  used  for  the  dipole 
inversion  (i.e.  we  only  bound  the  spatial  location).  For  scale  factors  we  use 

xtyp  =  (1, 1,  l)m,  (9typ,  (ptyp)  =  (7r/2,7r/2)  radians  and  dtyp  =  5nT  (34) 

To  determine  the  starting  model,  we  first  do  a  dipole  inversion  and  go  through  the  dipole 
identification  process.  For  each  ordnance  item,  the  identification  process  will  return  the 
dip  angle  9  and  azimuth  0  that  achieved  the  closest  fit  to  the  dipole.  These  can  be 
used  as  the  starting  orientation  with  the  dipole  location  used  as  the  starting  location, 
similarly  for  the  dc-shift.  Rather  than  doing  an  inversion  for  every  item  in  the  library 
we  only  consider  those  3-5  ordnance  items  with  the  best  dipole  fits. 
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The  Guthrie  Road,  Montana  dataset 

We  only  apply  the  dipole/octupole  identification  procedure  to  those  anomalies  flagged  as 
potential  UXO  by  the  remnant  magnetization  technique  (with  a  cutoff  of  50%  remnant 
magnetism).  The  quality  of  the  fits  obtained  for  the  optimum  ordnance  were  around 
the  same  as  the  dipole  fits;  sometimes  they  were  slightly  worse  and  sometimes  slightly 
better. 

Identification  results  for  the  dipole/octupole  method  for  76-mm  projectiles  are  given 
in  Figure  (12a)  and  for  81-mm  mortars  in  Figure  (12b).  For  76-mm  projectiles  the 
correct  identification  is  achieved  for  20  of  34  cases,  one  better  than  the  dipole  method. 
However,  for  the  81-mm  mortar  correct  identification  is  only  achieved  for  22  of  48  cases, 
three  less  than  the  dipole  method. 

Misidentification  results  for  false  alarms  are  given  in  Figure  (14b).  The  distribution 
of  items  is  very  similar  to  that  achieved  with  the  dipole  method. 

We  believe  that  there  are  three  reasons  why  the  dipole/octupole  identification  method 
does  not  improve  the  original  dipole  identification  method. 

1.  The  octupole  response  is  quite  small  and  does  not  provide  significant  additional 
constraints  on  the  the  geometry; 

2.  The  equivalent  spheroid  is  not  an  accurate  enough  approximation  of  an  ordnance 
item  for  the  octupole  contribution  to  have  a  positive  effect;  and 

3.  The  ordnance  items  are  slightly  remnantly  magnetized  which  changes  the  relation¬ 
ship  between  ordnance  orientation  and  the  direction  of  magnetization. 


Identification  using  the  dipole  and  quadrupole 

Thus  far,  the  identification  routines  we  have  developed  have  relied  solely  on  a  spheroid 
model  of  ordnance.  However,  this  assumption  is  overly  simplistic  since  most  ordnance 
items  do  not  have  that  degree  of  symmetry.  In  particular,  the  front  and  tail  ends  of 
ordnance  tend  to  be  different:  the  front  is  often  tapered  to  provide  better  aerodynamics 
and  target  penetration,  while  the  tail  end  is  often  fiat  (the  tails  of  mortars  don’t  matter 
because  they  are  usually  made  of  aluminum).  Therefore,  a  more  realistic  model  of 
ordnance  is  an  axially  symmetric  body  without  front-back  symmetry;  such  a  body  has 
a  multipole  expansion  with  both  a  dipole  term  and  a  non-zero  quadrupole  term  (recall 
that  the  quadrupole  was  zero  for  a  spheroid  due  to  symmetry). 

Our  intent  with  this  new  identification  method  is  to  try  to  use  the  information  in  the 
quadrupole  to  constrain  the  orientation  of  the  body  (</>,  9).  As  we  discovered  when  using 
the  dipole  method,  the  principal  source  of  ambiguity  in  identification  is  our  inability 
to  constrain  the  orientation  of  the  body.  In  particular,  when  the  object’s  dimension  is 
changed,  one  can  compensate  by  changing  the  object  orientation  and  produce  the  same 
or  a  very  similar  dipole  moment.  If  we  could  constrain  the  orientation  then  we  could 
better  constrain  the  ordnance  dimension. 

The  basis  of  the  idea  is  that  there  will  be  no  quadrupole  contributions  perpendicular 
to  the  axis  of  symmetry  but  there  will  be  contributions  along  the  symmetry  axis.  We 
therefore  construct  a  model  with  these  characteristics  and  use  an  inversion  procedure 
to  recover  a  model  with  the  best  orientation.  Concomitant  with  this  process  we  also 
recover  parameters  that  depend  on  the  shape  of  the  body,  as  we  now  explain. 
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The  dipole  moment  of  an  axially  symmetric  body  can  be  written  in  the  form, 


m  =  /x“1ArFAB0  (35) 

where  A  is  the  usual  Euler  rotation  tensor  and  F  is  a  3  x  3  diagonal  matrix  with 
(F2,  F2,  F$)  along  the  diagonal.  We  are  assuming  axial-symmetry  so  that  F\  =  F2. 
The  demagnetization  factors,  F),  given  in  the  above  equation,  are  VFl  times  those 
previously  given  in  Equation  2,  where  V  is  the  volume  of  the  body.  These  volume 
modified  demagnetization  factors  contain  all  the  information  on  the  body’s  dimensions 
that  we  are  attempting  to  recover. 

We  are  assuming  the  body  is  axially  symmetric  (about  the  z-axis),  therefore  the 
quadrupole  tensor  components  in  the  x-y  plane  are  zero  (see  Appendix  A  for  a  deriva¬ 
tion).  The  body  is  asymmetric  about  the  x-y  plane  so  that  the  quadrupole  tensor  has 

~  T 

non-zero  components  in  the  z-direction.  The  tensor  is  symmetric  Mf/  =  Mg  and  is 
therefore  defined  by  the  equation, 


0  0  mq  1 

0  0  mq  2 

mqi  mq  2  mq3 


(36) 


where  mq\,  mq 2  and  mq 3  are  dependent  on  the  shape  and  magnetization  of  the  body. 
The  above  equation  is  in  body-centered  coordinates.  The  quadrupole  tensor  may  be 
rotated  into  the  Earth’s  frame  of  reference  using  the  following  equation, 


M9  =  AMgAr  (37) 

The  magnetic  field  due  to  this  quadrupole  tensor  can  be  calculated  using  Equation  (A- 
11).  Ignoring  the  octupole  and  any  higher  order  moments,  full  specification  of  an  axially- 
symmetric  body  therefore  involves  the  following  11  parameters 


[x,  y,  z,  <f>,  9,  F2,  F3,  mgi ,  mq2,  mq3,  d\  (38) 

where  we  have  again  allowed  a  dc-shift,  d. 

To  invert  for  the  11  parameter  model  we  use  the  same  least-squares  formulation  and 
Interior-reflective  Newton  method  described  previously.  The  parameter  scalings  are  the 
same  as  those  given  in  Equation  (34),  with  the  additional  parameters  scaled  as 

F typ  =  (0.01,0.01)m3  and  mqtyp  =  (0.01,  .01,  .01)Am3  (39) 

Once  the  model  parameters  have  been  recovered  we  return  to  our  spheroid  model  of 
ordnance  and  find  the  object  dimensions  that  reproduce  the  recovered  volume  demag¬ 
netization  factors,  F2  and  F3.  This  is  straightforward  because  the  aspect  ratio,  e  solely 
determines  the  ratio  F2jF3,  while  the  diameter  and  aspect  ratio  together  determine  the 
magnitudes  of  F2  and  F3.  On  occasion  we  find  that  F2  converges  towards  zero  while  F3 
appears  to  converge  to  a  sensible  value.  In  that  case,  we  fix  the  aspect  ratio  at  a  typical 
value  3  <  e  <  5  and  just  solve  for  the  diameter. 

The  Guthrie  Road,  Montana  dataset 

We  applied  the  quadrupole  inversion  method  to  all  737  anomalies  at  Guthrie  Rd  that 
had  good  dipole  fits.  Perhaps  self-evidently  the  fits  and  correlation  coefficients  for 
the  dipole/quadrupole  model  were  better  than  for  the  dipole  model  alone.  Figure  15 
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Aspect  ratio  Ordnance  diameter  (mm) 


Figure  15:  Quadrupole  identification  results  at  Guthrie  Road;  (a)  Recovered  aspect- 
ratios  and  diameters;  and  (b)  Histogram  of  recovered  diameters. 


summarizes  the  results  of  the  quadrupole  identification.  There  appears  to  be  some 
identification  ability  but  there  are  a  significant  number  of  ordnance  items  where  the 
recovered  dimensions  are  inaccurate.  In  particular,  there  is  a  tendency  to  either  recover 
too  low  an  aspect  ratio  and  hence  too  large  a  diameter  or  vice-versa. 

There  are  two  conclusions  that  can  be  drawn  from  the  quadrupole  discrimination 
results.  Firstly,  the  quadrupole  may  be  unable  to  provide  the  required  constraints  on 
the  ordnance  orientation.  Alternatively,  there  may  be  enough  remnant  magnetization 
in  the  body  so  that  our  model  is  not  an  accurate  reflection  of  reality. 

With  these  issues  in  mind  we  trialed  one  last  method  for  identification  using  the 
quadrupole.  This  time  we  allow  the  dipole  moment  to  vary  independently  of  the 
quadrupole;  i.e.  in  place  of  F\  and  F2  in  Equation  (38)  we  solve  for  m  =  (7711,7712,7713). 
Thus  there  are  twelve  parameters  to  find.  Our  intent  is  to  determine  if  unmodelled 
quadrupole  components  bias  the  estimation  of  the  dipole  moment. 

We  applied  the  method  to  the  Guthrie  Road  data  and  then  used  the  recovered  dipole 
moment  to  achieve  ordnance  identification  in  the  usual  way.  We  found  that  the  method 
was  slightly  inferior  to  the  original  dipole  identification  routine.  In  particular  for  the 
81-mm  mortar  23  of  the  48  anomalies  were  identified  correctly  (compared  to  25  of  48), 
while  for  the  76-mm  projectile  only  16  of  34  cases  were  correct  (compared  to  20  of  34). 
Thus,  we  conclude  that  unmodelled  quadrupole  components  do  not  negatively  affect  our 
ability  to  identify  using  the  dipole  moment. 
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7  Error  Analysis 


The  dipole  identification  results  for  Guthrie  Road  indicate  that  around  55%  of  the 
time  we  correctly  identify  the  ordnance  type.  The  45%  misidentifications  must  be  due 
to  either  an  incorrect  model,  remnant  magnetism  or  uncertainty  in  the  recovered  dipole 
moments.  In  this  section  we  analyze  the  uncertainty  in  recovered  dipole  moment  using 
linearized  error  analysis. 

Estimating  the  covariance  matrix  of  the  model  parameters 

To  derive  the  covariance  matrix  of  the  model  parameters  we  follow  the  development 
given  in  Bard  (1974).  At  the  minimum  of  the  objective  function  the  gradient  is  zero 

V0(m,d)  =  0  (40) 

where  V  is  with  respect  to  the  model  parameters.  If  we  change  the  data  to  d  +  dd  we 
will  change  the  objective  function  and  hence  the  position  of  the  minimum  will  vary  to 
m  +  dm.  The  gradient  of  the  new  objective  function  at  the  new  minimum  will  be  zero. 

V<£(m  +  dm,  d  +  dd)  =  0  (41) 

If  we  expand  this  gradient  in  a  Taylor  series  and  only  keep  first  order  terms  we  find 

A 

V0(m  +  dm,  d  +  dd)  =  V<?i(m,  d)  +  V2^(m,  d)dm  +  — V(£(m,  d)dd  —  0  (42) 

o  d 

Recognizing  that  the  V2  term  is  the  Hessian  and  rearranging  we  find 

A 

dm  =  H-1  V0(m,  d)dd  (43) 

ad 

The  covariance  matrix  of  the  model  parameters  is 

Vm  =  E[dmdmr]  (44) 

where  E[-]  is  the  expectation  operator.  Substituting  Equation  (43)  into  (44),  recognizing 
that  the  Hessian  and  gradient  terms  are  constants  (under  expectation)  and  rearranging 
terms  we  find 


A  A 

Vm  =  H-1— V^(m,d)rE[ddddT]— V<^(m,d)H-1  (45) 

The  covariance  matrix  of  the  data  is  =  E[ddddT].  Assuming  the  data-weighting 
matrix  in  our  formulation  of  the  inverse  problem  is  the  same  as  the  inverse  covariance 
matrix  of  the  data,  the  Gauss-Newton  approximation  to  the  Hessian  is 

H  =  JTVd1J  (46) 

while 


— V0(m,d)=Vd-1J 


(47) 
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Substitution  of  Equations  (46)  and  (47)  into  Equation(45)  and  simplification  reveals 
that  the  covariance  matrix  is  the  inverse  of  the  Gauss-Newton  approximation  to  the 
Hessian, 


Vm  =  (JT VJ)-1  -  H"1  (48) 

The  above  equation  can  be  used  to  estimate  uncertainties  in  each  of  the  model  param¬ 
eters  and  can  also  be  used  to  define  confidence  ellipses  of  multiple  parameters.  The 
uncertainties  obtained  through  this  equation  are  linearized  estimates  of  the  true  un¬ 
certainties.  Before  considering  them  further  we  investigate  how  closely  these  linearized 
estimates  would  match  fully  non-linear  estimates.  If  the  minimum  of  the  objective 
function  is  m  then  linearized  estimates  are  equivalent  to  approximating  the  objective 
function  with 


0(m)  —  ^>(m)  «  ^(m  -  m)TH(m  -  m)  (49) 

If  the  above  equation  is  a  good  approximation  to  the  objective  function  about  the 
minimum  then  the  linearized  uncertainties  will  be  accurate  estimates  of  the  true  uncer¬ 
tainties.  We  have  compared  the  objective  functions  of  many  anomalies  and  usually  find 
a  situation  such  as  that  illustrated  in  Figure  (16).  The  linear  and  non-linear  estimates 
are  very  similar  for  parameters  related  to  the  dipole  moment,  and  less  so  for  parameters 
related  to  location  (especially  the  depth).  From  the  perspective  of  ordnance  identifica¬ 
tion  it  is  the  uncertainties  in  the  dipole  moment  that  are  important,  so  we  can  conclude 
that  the  linearized  estimates  should  be  accurate  for  our  purposes.  We  note  in  passing 
that  including  the  second  order  terms  in  the  Hessian  has  very  little  impact  on  the  accu¬ 
racy  of  the  linear  approximation,  even  if  the  residuals  at  the  solution  are  large.  This  is 
readily  apparent  for  the  dipole  components  in  cartesian  coordinates  m  =  (mx,my,mz) 
because  the  second  derivatives  with  respect  to  these  parameters  are  zero  (Equation  9). 

Confidence  ellipses  of  multiple  parameters 

We  are  principally  interested  in  the  joint-uncertainty  in  the  dipole  magnitude  and  an¬ 
gle  relative  to  the  Earth’s  field  because  these  two  parameters  control  our  decisions  on 
discrimination  and  identification.  A  good  method  for  visualizing  the  joint-uncertainties 
is  to  calculate  confidence  ellipses;  these  are  contours  of  equal  confidence  in  the  model 
parameters.  They  are  usually  defined  to  encompass  a  certain  probability  that  the  model 
lies  within  the  defined  region.  For  instance,  a  one-sigma  (or  68.3%)  confidence  ellipse 
would  delineate  a  region  that  has  a  68.3%  chance  of  containing  the  true  model  parame¬ 
ters  (under  the  particular  statistical  assumptions  used  to  define  the  objective  function). 

To  define  the  confidence  ellipse  for  r  parameters  requires  that  the  uncertainties  in 
the  other  n  —  r  parameters  (assuming  there  are  n  parameters  in  total)  are  projected 
onto  the  r-dimensional  hyperplane  under  consideration.  In  turns  out  (Bard,  1974)  that 
this  projection  can  be  implicitly  achieved  by  calculating  the  covariance  matrix  through 
inversion  of  the  appropriate  r  x  r  subsection  of  the  Hessian  matrix.  For  example,  if  we 
are  interested  in  the  joint  confidence  ellipse  for  parameters  2  and  5  we  would  calculate 

-l 

(50) 

The  principal  axes  of  the  ellipse  are  then  obtained  by  eigenvalue  decomposition  of  Vm, 
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Figure  16:  Objective  function  (solid  line)  and  Gauss-Newton  approximation  to  the  Hes¬ 
sian  (dashed  line)  about  the  inversion  solution  for  Target  575,  Fort  Ord. 


Vm  =  USUr  (51) 

where  S  is  a  diagonal  matrix  of  eigenvalues  and  U  is  an  orthonormal  matrix  of  eigen¬ 
vectors;  each  eigenvector  is  a  principal  axis  of  the  ellipse.  The  size  of  each  principal  axis 
is  given  by  the  reciprocal  of  the  corresponding  eigenvalue. 

The  final  step  is  to  find  the  ellipse  that  lies  on  the  boundary  of  a  specified  confidence 
region  (e.g.  68.3%  or  the  one-sigma  error).  If  the  sampling  distribution  is  normal,  unbi¬ 
ased  and  we  assume  that  the  covariance  matrix  Vm  is  known,  then  (m— m)TV“1(m— m) 
is  x2  distributed  with  r  degrees  of  freedom  (Bard,  1974).  The  boundary  of  the  ellipse 
corresponding  to  a  particular  confidence  level  is  found  by  calculating  the  value  of  the  y2 
distribution  with  r  degrees  of  freedom  that  corresponds  to  that  confidence.  For  example, 
with  2  degrees  of  freedom,  the  68.3%  confidence  limit  has  x2  =  2.298.  The  contour  of 
(m  —  m  )Tv-1  (m  —  m)  with  the  value  2.298  corresponds  to  the  boundary  of  the  68.3% 
confidence  region. 

Alternatively,  if  the  covariance  matrix  is  not  known  and  we  assume  that  all  observa¬ 
tions  are  independent  (so  that  V<f  =  <r2I)  then  the  confidence  level  is  obtained  through 
the  FT)Tl-r  distribution  (Bard,  1974).  In  that  case  the  value  of  x2  calculated  above  would 
be  replaced  with 


ra2Fr)„_r(  7) 


(52) 
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where  7  is  the  desired  confidence  limit  (7  =  0.683  for  the  68.3%  confidence  level)  and  cr2 
is  the  sample  variance.  If  the  optimum  model  returned  by  the  minimization  routine  is  m 
and  there  are  K  model  parameters  then  the  sample  variance  is  obtained  by  calculating 

<’2  =  iv^£[Fi(,h)“<^]2  (53) 

1=1 

Finally,  if  the  assumption  of  Gaussian  statistics  does  not  hold,  conservative  confidence 
limits  can  be  calculated  by  recourse  to  the  Bienayme-Chebyshev  inequality  (Bard,  1974). 
For  the  7  confidence  limit  this  replaces  the  value  of  y2  above  with 

r 
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The  Guthrie  Road,  Montana  dataset 

The  distribution  of  residuals  for  the  660  anomalies  at  Guthrie  Road  with  correlation 
coefficients  greater  than  0.8  are  shown  in  Figure  (17).  The  curves  have  been  normalized 
so  that  they  have  an  integral  of  one;  thus  we  are  approximating  the  posterior  probability 
density  function  (posterior  PDF)  of  the  residuals.  The  best  fitting  Gaussian  distribution 
(with  a  =  0.99  nT)  provides  a  poor  fit  to  the  residuals.  In  particular,  large  residuals  are 
much  more  likely  than  predicted  by  a  normal  distribution.  In  Figure  (17)  we  also  show 
the  best  fitting  Ekblom  distribution  (Ekblom,  1973)  which  has  a  probability  density 
function  equal  to 


Pr(x)  =  cexp 


(x2  +  e2)^/2 
aP 


(55) 


where  c  is  a  constant  that  ensures  the  PDF  integrates  to  unity.  This  is  a  perturbed 
version  of  an  Lp-norm  and  we  found  best  fitting  values  of  a  =  1.16  x  10-4  nT,  e  — 
0.609  nT  and  p  =  0.224.  The  correspondence  between  the  observed  distribution  of 
residuals  and  the  Ekblom  PDF  is  excellent.  We  will  consider  the  good  fit  of  the  Ekblom 
PDF  further  in  the  discussion  section  of  this  report. 

As  the  residuals  are  clearly  non-Gaussian,  we  must  be  cautious  when  interpreting 
confidence  ellipses  calculated  under  the  normality  assumption.  The  shape  and  relative 
size  of  the  ellipses  are  likely  to  be  correct  but  their  absolute  values  are  probably  incorrect 
(i.e.  the  68.3%  confidence  ellipse  won’t  correspond  to  a  confidence  of  68.3%).  However, 
the  absolute  confidence  ellipses  should  provide  a  rough  estimate  of  the  uncertainties 
in  the  parameters.  They  should  also  provide  a  reliable  means  of  comparison  of  the 
uncertainties  in  recovered  parameters  for  different  anomalies. 

We  calculated  confidence  ellipses  for  the  dipole  magnitude  and  angle  relative  to  the 
Earth’s  field,  for  the  anomalies  due  to  76-mm  projectiles  and  81-mm  mortars  at  Guthrie 
Road  (Figure  18).  The  data- weighting  matrix  was  the  identity  matrix,  with  a  sample 
standard  deviation  calculated  from  the  residuals  by  Equation  (53).  The  figures  show 
the  one-sigma  (68.3%)  and  two-sigma  (95.4%)  confidence  ellipses  about  each  recovered 
dipole  moment  (assuming  Gaussian  statistics).  For  the  76-mm  projectiles  all  except  one 
anomaly  have  very  small  uncertainties.  The  one  outlier  was  incorrectly  identified  as  a 
105-mm  projectile.  For  the  81-mm  mortars  all  anomalies  have  very  small  uncertainties. 
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Figure  17:  Distribution  of  residuals  for  the  660  anomalies  at  Guthrie  Rd  with  correlation 
coefficient  greater  than  0.8. 

We  also  calculated  confidence  ellipses  using  the  very  conservative  Bienayme-Chebyshev 
inequality  (Figure  19).  The  confidence  ellipses  axe  significantly  larger  than  those  re¬ 
turned  under  the  normal  assumption.  If  the  real  confidence  ellipses  were  this  large,  then 
uncertainty  in  the  recovered  moment  would  be  a  contributing  factor  in  misidentifications 
(especially  for  76-mm  projectiles).  However,  the  real  confidence  ellipses  would  have  a 
size  intermediate  between  the  Gaussian  and  Bienayme-Chebyshev  estimates.  This  would 
imply  that  misidentification  is  sometimes,  though  probably  not  that  often,  caused  by 
uncertainty  in  the  recovered  moment.  Misidentifications  are  more  likely  to  occur  due  to 
modelling  deficiencies  and/or  remnant  magnetization. 

We  also  calculated  confidence  ellipses  (under  the  normal  assumption)  for  the  anoma¬ 
lies  due  to  metallic  debris  and  geology  (Figure  20).  For  the  metallic  debris  many  of  the 
uncertainties  are  quite  low  but  there  are  a  few  items  with  much  higher  uncertainties. 
For  the  anomalies  due  to  geology  the  uncertainties  are  generally  much  higher.  In  many 
cases  the  recovered  angle  is  around  60°  but  the  uncertainty  is  high.  This  is  a  critical  area 
for  making  discrimination  decisions  so  we  can  conclude  that,  for  geological  anomalies, 
uncertainty  in  dipole  angle  can  have  a  significant  impact  on  discrimination. 
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Figure  20:  One  sigma  (68.3%)  and  two  sigma  (95.4%)  confidence  ellipses  about  the 
recovered  dipole  moments  for  (a)  metallic  debris;  and  (b)  geological  anomalies. 
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8  Discussion 


In  this  report  we  have  described  a  successful  inversion  method  for  recovering  the 
position  and  dipole  moment  from  a  total-field  magnetic  anomaly.  The  inversion  uses 
a  least-squares  formulation  and  an  interior-reflective  Newton  method  to  minimize  the 
objective  function.  A  starting  model  is  defined  using  an  analytical  method  based  on 
approximating  the  vertical  component  of  the  magnetic  field.  The  inversion  method  is 
very  fast  and  reliable;  our  Matlab  implementation  requires  only  a  few  seconds  to  recover 
the  parameters  for  each  anomaly. 

The  recovered  dipole  moments  can  be  used  to  make  discrimination  decisions  by 
recourse  to  the  phenomena  of  shock  demagnetization.  Complete  shock  demagnetization 
may  not  occur  for  all  ordnance  items  all  the  time.  Additionally,  non-ordnance,  even 
though  they  may  be  remnantly  magnetized,  might  be  orientated  so  that  their  dipole 
moments  lie  close  to  the  Earth’s  field.  This  means  that  it  is  not  possible  to  make  a 
definitive  ordnance/non-ordnance  decision.  Rather,  we  prefer  to  rank  anomalies,  with 
items  most  likely  to  be  UXO  higher  up  the  list.  We  investigated  three  such  ranking 
procedures;  (i)  dipole  angle,  (ii)  dipole  angle  and  magnitude  and  (iii)  percentage  of 
remnant  magnetism,  and  found  that  the  remnant  magnetization  ranking  was  the  best. 
With  this  method  all  ordnance  would  have  been  recovered  after  only  half  of  the  anomalies 
had  been  excavated.  This  included  85  anomalies  with  poor  dipole  fits  that  only  yielded 
one  UXO.  The  criterion  used  to  select  potential  UXO  may  need  to  be  revised  to  minimize 
the  number  of  anomalies  where  we  are  unable  to  make  reliable  inferences  about  the 
bodies  origin. 

Linearized  error  analysis  indicates  that  the  angle  between  the  recovered  dipole  mo¬ 
ment  and  the  Earth’s  field  for  ordnance  items  is  relatively  well  constrained.  Therefore, 
uncertainty  in  the  recovered  dipole  moment  is  unlikely  to  strongly  influence  the  dis¬ 
crimination  ranking.  For  items  identified  as  geology,  the  confidence  ellipses  were  much 
larger.  Consequently,  the  position  of  these  items  in  the  ranking  is  strongly  influenced 
by  uncertainty. 

The  discrimination  method  was  very  successful  at  Guthrie  Road  but  could  not  be 
reliably  applied  to  the  seeded  sites  at  Fort  Ord.  This  is  a  consequence  of  remnant  mag¬ 
netization  and  emphasizes  the  importance  of  shock  demagnetization  to  the  viability  of 
the  method.  Additional  work  needs  to  be  undertaken  to  determine  the  circumstances 
where  shock  demagnetization  occurs.  For  instance,  the  recovered  dipole  moments  of 
several  of  the  81-mm  mortars  at  Guthrie  Road  had  relatively  large  dipole  angles,  in¬ 
dicating  the  presence  of  some  remnant  magnetism.  Perhaps  the  lower  impact  velocity 
of  mortars  compared  to  artillery  projectiles  causes  only  partial  shock  demagnetization; 
this  needs  to  be  investigated  further. 

Identification  of  the  ordnance  type  using  the  recovered  dipole  moment  is  difficult  due 
to  non-uniqueness.  The  dipole  moment  is  the  product  of  volume  with  the  magnetization 
vector.  When  the  volume  changes,  self-demagnetization  as  the  orientation  of  the  body 
is  varied  can  change  the  magnetization  vector  and  thus  produce  an  identical  or  similar 
dipole  moment.  However,  recognizing  that  there  are  only  a  finite  number  of  ordnance 
types  in  an  given  area,  makes  identification  using  the  recovered  dipole  moment  feasible. 
At  Guthrie  Road  correct  identification  could  be  achieved  about  55%  of  the  time. 

One  of  the  key  factors  that  influences  successful  identification  is  obtaining  the  correct 
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trajectory  of  moment  magnitude  and  angle  for  a  given  ordnance  item.  We  use  a  spheroid 
approximation  to  calculate  this  trajectory.  Modelling  of  shapes  that  more  closely  ap¬ 
proximate  real  ordnance  should  be  undertaken  in  order  to  determine  the  reliability  of 
the  spheroid  approximation. 

A  significant  cause  of  misidentification  is  the  inability  of  the  dipole  moment  to  con¬ 
strain  the  orientation  of  the  causative  body.  In  an  attempt  to  overcome  this  problem  we 
developed  inversion  routines  that  recover  information  on  the  quadrupole  and  octupole 
components  of  the  body.  The  dipole/octupole  method  we  developed  uses  a  spheroid 
approximation  of  ordnance  to  estimate  both  the  dipole  and  octupole  moments.  No  sig¬ 
nificant  improvement  in  identification  over  the  dipole  method  was  found.  This  could 
be  due  to  one  or  more  of  the  following:  (i)  remnant  magnetization;  (ii)  rapid  falloff  of 
the  octupole  component  with  distance  so  that  it  is  unable  to  constrain  the  orientation; 
and  (iii)  the  spheroid  provides  a  poor  approximation  to  the  octupole  moment  of  real 
ordnance. 

The  quadrupole  method  relaxed  the  requirement  of  top-bottom  symmetry  enforced 
by  a  spheroid  approximation  and  attempted  to  recover  the  dipole  and  quadrupole  mo¬ 
ments  of  an  arbitrary  axially  symmetric  body.  The  fits  to  the  data  were  improved 
slightly  but  no  advantages  for  identification  were  realized.  The  cause  is  likely  due  to  one 
or  more  of  the  following:  (i)  remnant  magnetization;  (ii)  the  inability  of  the  quadrupole 
component  to  constrain  the  orientation;  and  (iii)  inadequacy  of  the  model. 

The  failure  of  the  quadrupole  and  octupole  methods  to  improve  identification  may 
imply  that  analysis  of  magnetic  data  alone  does  not  provide  enough  information  to 
constrain  the  objects  orientation.  One  possible  method  to  provide  this  constraint  is  by 
joint  inversion  of  magnetic  and  electromagnetic  data;  either  single  or  multiple  channel 
time  domain  EM  or  frequency  domain  EM.  Most  EM  sensors  illuminate  the  target  from 
multiple  directions  and  are  consequently  more  sensitive  to  orientation  (in  contrast  the 
direction  of  the  Earth’s  field  is  fixed).  In  return,  the  magnetics  may  provide  valuable 
depth  constraints  that  are  difficult  to  obtain  with  EM  alone  (Pasion  and  Oldenburg, 
2001). 

Linearized  error  analysis  was  used  to  estimate  individual  parameter  uncertainties  and 
joint  confidence  ellipses.  The  linear  approximation  provides  an  accurate  model  of  the 
objective  function  about  the  minimum.  However,  the  residuals  do  not  follow  a  Gaussian 
distribution  so  that  reliable  confidence  ellipses  could  not  be  defined.  Additionally,  the 
upper  bound  provided  by  the  Bienayme-Chebyshev  inequality  is  very  conservative  and 
probably  returns  confidence  ellipses  that  are  much  larger  than  required.  A  recommen¬ 
dation  for  future  work  is  to  attempt  to  calculate  confidence  ellipses  using  non-Gaussian 
statistics. 

The  residuals  were  found  to  follow  an  Ekblom  distribution,  yet  we  minimized  a 
measure  derived  from  Gaussian  statistics.  A  reformulation  of  the  problem  using  an 
Ekblom  distribution  as  a  measure  of  data  misfit  is  therefore  recommended.  The  resulting 
objective  function  would  need  to  be  solved  by  a  different  optimization  algorithm  than 
the  one  used  here.  Algorithms  for  the  solution  of  this  optimization  problem  are  well 
known  and  include  iteratively  reweighted  least  squares  and  Newton’s  method  (Watson, 
1980;  Ekblom,  1987).  Implementation  of  the  new  optimization  method  would  enable  a 
thorough  investigation  to  be  made  on  the  influence  the  choice  of  the  measure  of  data 
misfit  has  on  the  recovered  model.  This  is  an  issue  often  glossed  over  when  solving 
inverse  problems  and  this  report  was  no  exception. 
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9  Conclusions  and  Recommendations 


A  dipole  inversion  method  was  developed  and  successfully  applied  to  the  discrim¬ 
ination  and  identification  of  magnetic  anomalies  at  Guthrie  Road  Montana.  The  dis¬ 
crimination  method  had  the  potential  to  reduce  the  number  of  excavations  by  half,  yet 
still  yielded  all  the  ordnance.  Identification  of  the  causes  of  the  anomalies  is  difficult 
using  the  recovered  dipole  moment  due  to  non-uniqueness.  However,  by  recognizing  that 
there  are  a  finite  number  of  ordnance  types,  correct  identification  could  be  achieved  in 
about  55%  of  cases.  Extension  of  the  inversion  method  to  recover  either  quadrupole  of 
octupole  moments  did  not  improve  our  ability  to  discriminate  or  identify. 

We  recommend  the  following  areas  as  directions  for  future  research. 

1.  Shock  demagnetization:  The  extent  to  which  this  phenomena  occurs  should  be 
investigated  further  because  shock  demagnetization  is  central  to  the  success  of 
discrimination  using  magnetic  data. 

2.  Magnetic  modelling  of  ordnance:  The  trajectory  of  the  modelled  dipole  moment 
magnitude  and  angle  determines  our  identification  decisions.  Therefore,  the  agree¬ 
ment  between  the  moments  returned  by  a  spheroid  approximation  and  a  more 
realistic  ordnance  shape  should  be  investigated  further. 

3.  Joint  inversion  of  magnetics  and  electromagnetics:  Identification  difficulties  in 
magnetics  are  caused  by  an  inability  to  constrain  the  orientation  of  the  causative 
body.  The  ability  of  time  or  frequency  domain  electromagnetics  to  constrain  the 
orientation  needs  to  be  investigated.  Conversely,  magnetics  may  provide  valuable 
depth  information  that  can  be  exploited  by  the  EM  inversion. 

4.  Error  analysis  using  non-Gaussian  statistics:  The  residuals  do  not  follow  a  Gaus¬ 
sian  distribution  so  uncertainty  estimates  based  on  a  more  realistic  statistical 
distribution  should  be  investigated. 

5.  Robust  measures  of  data  misfit:  The  problem  should  be  reformulated  using  a  more 
realistic  measure  of  data  misfit  (the  Ekblom  norm),  and  methods  implemented  to 
solve  the  modified  inversion  problem. 
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Appendix  A:  Multipole  Expansion 


In  this  appendix  we  give  a  brief  derivation  of  the  multipole  expansion  of  a  magnetic 
body,  and  in  particular  derive  the  equations  relevant  to  a  spheroid.  Our  derivations 
closely  follow  those  given  in  McFee  (1989). 

The  scalar  potential,  <f>( r)  of  a  body  with  magnetization  M  is  given  by 

where  r'  is  the  distance  from  the  source  to  the  observation  point  (Figure  1). 
the  gradient  and  using  the  divergence  theorem  we  find 

*V  =  Ws?MdS-Jv?V-MdV)  <A-2> 

If  the  induced  field  inside  the  body  is  uniform  and  parallel  (which  it  is  for  a  spheroid), 
which  eliminates  the  second  term  in  the  above  equation  because  then  V  •  M  =  0.  The 
multipole  method  proceeds  by  expanding  ^  as  a  Taylor  series  about  the  origin, 


(A-l) 
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where  repeated  indices  imply  that  summation  convention  is  to  be  used.  The  scalar 
potential  from  a  2"  pole  with  moment  is  given  by  the  expression 

(A-4) 

from  which  it  follows  that  the  first  four  moments  are 
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(A-7) 

(A-8) 


In  order  the  moments  are  the  scalar  monopole  the  vector  dipole  m^\  the  rank  2 
tensor  quadrupole  m!2\  and  the  rank  3  tensor  octupole  Magnetic  monopoles  do 


A1 


not  exist  which  means  that  the  monopole  term  is  zero.  Expansion  of  the  dipole  term 
reveals  that  the  dipole  moment  is  the  product  of  the  magnetization  with  the  volume, 

m*1)  =  ME  (A-9) 

Due  to  the  symmetry  of  the  spheroid  the  quadrupole  term  is  zero.  After  the  dipole,  the 
next  non-zero  term  is  the  octupole;  specification  of  the  octupole  components  is  tedious 
so  we  refer  the  interested  reader  to  McFee  (1989). 

The  magnetic  field  (in  the  i-th  direction)  due  to  the  dipole  component  is  given  by 
the  expression 


Bl1]  =  ^3  (^tX  '  m(1)]X«  -  m‘)  (A‘10) 

Although  zero  for  the  spheroid  we  will  use  the  quadrupole  term  later,  its  field  is, 

B?]  =  ^  (~ximjf  -  xAmi?  +  mjf)  +  5r-2xiXjXkmff')  (A-ll) 
and  for  the  octupole  the  field  is, 


•B,(3)  =  "5  (3n»y]  -  15r  2  +  XjXkm{^  +  35r  ixlxjxkxlmf^  (A-12) 

It  remains  to  find  the  magnetization  vector  for  a  spheriod.  A  solution  of  the  boundary 
value  problem  (Stratton,  1941)  shows  that  the  demagnetization  factors  are 


p  — _ _ * _ 

1  +  ai(/ir  -  l)/2 

where  /ir  is  the  relative  permeability  (p  =  pT/i0), 
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for  a  prolate  spheroid  (e  >  1)  and 


arctan(7Tif )  “  */2 
\/l  —  e2 


for  an  oblate  spheroid  (e  <  1).  For  the  special  case  of  a  sphere  (e  =  1)  then 


ai  =  «2 


(A-13) 
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The  demagnetization  factors  together  with  the  Earth’s  field  B0  (in  spheroid  centered 
coordinates),  determine  the  strength  of  the  induced  magnetization, 

M  =  ft-'FBo  (A- 19) 

where  we  have  written  F  as  a  3  x  3  diagonal  matrix  with  (F2,  F2,  F3)  along  the  diagonal 
(because  F\  —  F2). 


Multipole  expansion  for  an  axi-symmetric  body 

In  going  from  the  general  equation  (A-2)  for  a  magnetized  body  to  that  for  a  spheroid, 
we  used  the  fact  that  the  magnetization  inside  the  body  was  constant  so  that  V  •  M  =  0. 
For  a  general  axi-symmetric  body  this  is  no  longer  true  and  there  is  an  extra  contribution 
to  each  of  the  moments.  For  the  dipole, 
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and  using  the  identity 
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we  find 
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Thus  the  dipole  moment  is  just  the  integration  of  the  magnetization  over  the  volume 
of  the  body.  Most  ordnance  items  are  either  cylinder  like  with  a  pointed  front  (e.g. 
most  artillery  projectiles)  or  close  to  a  spheroid  shape  (e.g.  mortars,  the  fins  can  be 
ignored  as  they  are  usually  aluminum).  For  an  infinite  cylinder  the  field  lines  are  con¬ 
stant  within  the  body.  Therefore,  for  a  projectile  any  changes  in  magnetization  will  be 
concentrated  at  the  ends  of  the  body,  and  will  generally  make  only  a  small  contribution 
to  the  multipole  moments.  For  a  mortar,  the  shape  is  close  to  second  order  (Stratton, 
1941)  and  hence  the  magnetization  direction  will  be  approximately  constant.  Thus, 
changes  in  magnetization  direction  should  not  make  a  significant  contribution  to  the 
observed  magnetic  anomaly. 

Using  an  analogous  development  to  the  dipole,  the  quadrupole  moment  can  be  shown 
to  be  equal  to 


m\f  =  J^  M  ■  V(xiXj)dV 

=  J  ( MiXj  +  MjXi)dV  (A-23) 

Assuming  that  changes  in  magnetization  can  be  ignored  implies  that 

=  Mi  dx 3  [  XjdS  +  M7-  [  dx 3  f  XidS  (A-24) 
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where  L  is  the  length  of  the  body  parallel  to  the  symmetry  axis,  z  —  x 3  and  S(x 3)  is  a 
slice  of  the  body  perpendicular  to  the  axis  of  symmetry.  Due  to  symmetry  each  of  the 
integrals  over  S(x 3)  are  zero  for  i  or  j  =  1  or  2,  from  which  it  follows  that 

=  0  for  any  i,  j  —  1  or  2  (A-25) 

(2)  (2) 

The  quadrupole  moment  is  symmetric  rrv^  =  so  that  the  full  moment  tensor  can 
then  be  written  as  (where  we  have  dropped  the  (2)  subscript), 

0  0  mi3 

M  =  0  0  m23  (A-26) 
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Appendix  B:  Constraints  on  Dipole  Orientation 


The  angle  the  induced  dipole  moment  in  a  spheroid  makes  with  the  Earth’s  field 
cannot  exceed  a  certain  value.  In  this  appendix  we  derive  the  equations  which  specify 
this  upper  bound.  We  assume  the  Earth’s  field  is  at  an  angle  9  to  the  spheroid  symmetry 
axis  and  align  the  z-axis  of  our  coordinate  system  with  the  symmetry  axis.  We  can  then 
choose  the  x-axis  so  that  it  is  perpendicular  to  the  Earth’s  field.  In  this  coordinate 
system 


bo  =  bo(0,  sin  9 ,  cos  0)T  (B-l) 

The  spheroid  diameter  a  and  aspect  ratio  e  determine  the  demagnetization  factors,  F2 
and  T3,  so  that  the  induced  dipole  moment  m  is, 

m  —  V  ba( 0,  F2  sin  0,  F3  cos  0)T  (B-2) 

where  V  is  the  volume  of  the  spheroid.  The  angle  ip  between  the  Earth’s  field  and  the 
dipole  moment  is 


cos  ip  — 


m 


llbolllMI 

which,  on  substitution  of  Equations  (B-l)  and  (B-2),  becomes 

F2  sin2  9  +  F3  cos2  9 


cos  ip  = 


\J F$  sin2  9  +  F%  cos2  9 


(B-3) 
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After  differentiating  this  equation  by  9,  setting  the  result  to  zero  and  extensive 
algebraic  manipulation  we  find 


9  —  arctan 


(  I  Fj  +  F3Fj  -  2FWj\ 

VV  Fi  +  F2Fi  -  2F3F| J 


(B-5) 


This  orientation  of  the  spheroid  symmetry  axis  relative  to  the  Earth’s  field  makes  the 
angle  between  the  induced  dipole  and  the  Earth’s  field  a  maximum.  This  maxirrnmi 
angle  is  obtained  by  substituting  Equation  (B-5)  into  Equation  (B-4). 


B1 


Appendix  C:  Ambiguity  in  the  Dipole  Solution 


Once  the  sensor  distance  exceeds  a  few  body  lengths  distance  away  from  the  spheroid, 
the  field  is  essentially  dipolar.  In  such  cases,  inversion  for  spheroid  dimensions  has  an 
inherent  ambiguity:  many  different  spheroid  diameters  and  aspect-ratios  can  reproduce 
the  observed  dipole.  To  characterize  this  ambiguity  we  chose  our  coordinate  system  so 
that  the  z-axis  is  aligned  with  the  Earth’s  field, 


Bo  =  (0,0,Bo)t  (C-l) 

We  want  to  determine  if  the  dipole  m  from  a  particular  spheroid  (diameter  a  and 
aspect-ratio  e)  can  match  an  observed  dipole  moment  m.  Now  the  induced  dipole  for 
the  spheroid  is  given  by  the  following  equation 


m  =  VAtFABq 

Expanding  this  equation  we  find 
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(F2  —  F3 )  cos  9  sin  9  cos  ip 
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Now  we  want  to  find  a,  e,  9  and  ip  so  that  m  =  m.  We  can  find  an  expression  for  the 
spheroid  diameter  by  equating  the  norms  of  the  two  dipole  moments,  resulting  in 
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This  implies  that  the  spheroid  has  the  same  azimuth  as  the  dipole  moment. 

Determination  of  9  is  a  little  more  difficult.  We  start  by  using  the  second  row  of 
Equation  (C-3)  and  substitute  in  Equation  (C-4)  for  the  diameter, 


m2  _  (i*2  -  F3)  cos  9 sing  (n  ^ 

cos^IHI  “  yj F$  cos2 9  +  sin2  9 

Note  that  if  cos  ip  =  0  we  would  use  the  first  row  of  Equation  (C-3)  to  derive  an  equation 
in  terms  of  m\  and  sin  ip  .  Let  the  LHS  of  Equation  (C-7)  be  w,  expand  out  the  square- 
root  on  the  denominator  of  the  RHS,  and  use  the  substitution 


sin2  9  =  1  —  cos2  9 


(C-8) 
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After  some  algebraic  manipulation  Equation  (C-7)  reduces  to 


(F2  -  Fif  cos4  9  [a ?{Fl  -  F|)  -  {F2  -  F3)2]  cos2  9  +  oj2F2  (C-9) 

This  is  a  quadratic  equation  in  cos2  9  which  has  a  real  solution  if 

A  =  S2  -  4oc  >  0  (C-10) 

where 


a  =  (F2-  F3)2  (C-ll) 

b  =  -  Fi)  -  (F2  -  F3)2  (C-12) 

c  =  u>2F$  (C-13) 


The  solution  in  terms  of  6  is  then 


9  =  ±  arccos 


(C-14) 


The  ±  terms  in  the  above  equation  emphasize  that  when  A  >  0  there  are  actually  two 
possible  solutions. 

Determination  of  which  spheroids  can  reproduce  the  observed  dipole  proceeds  as 
follows.  Firstly,  we  chose  an  aspect-ratio  e  and  calculate  the  demagnetization  factors  in 
the  usual  way.  We  then  substitute  Equations  (C-ll)  to  (C-13)  into  Equation  (C-10)  to 
determine  if  a  solution  can  be  found  for  this  aspect-ratio.  If  it  can,  we  calculate  9  by 
Equation  (C-14),  and  finally  we  calculate  the  spheroid  diameters  (usually  there  will  be 
two)  by  Equation  (C-4). 
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